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‘^higher  current  density  obtained  in  the  second  data  set  was  due  to  the 
occurrence  of  the  co-precipitation  reaction  between  cadmium  ions  and  hy¬ 
droxy  ions  which  were  generated  by  the  electrochemical  reaction. 


Two  models,  which  described  the  transport  process  and  reactions  under 
two  different  experimental  conditions  were  developed.  The  model  for  a 
cadmium-ion-free  electrolyte  was  used  to  Identify  the  rate  of  the 
heterogeneous  electrode  reaction  by  matching  the  model  prediction  with 
the  experimental  results.  The  heterogeneous  reaction  rate  constants  so 
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CHAPTER  1 

INTRODUCTION 

Nickel-cadnium  storage  catteries  have  been  in  use  for 
more  than  fifty  years;  since  Norld  <Var  II,  new  design  princi¬ 
ples  have  been  introduced  which  have  greatly  extended  their 
usefulness.  Many  theories  and  investigations  have  been  put 
forward  concerning  the  reaction  mechanisms  of  both  cadmium  and 
nickel  electrodes.  Because  the  activity  of  the  nickel  electrode 
has  been  more  fully  characterized  than  that  of  the  cadmium 
electrode,  we  will  focus  our  attention  on  the  behavior  of  the 
cadmium  electrode. 

There  are  a  number  of  methods  used  for  the  fabrication 

of  cadmium  electrodes.  These  methods  are  usually  classified 

into  three  categories:  (a)  pressed  powder  or  pasted  type’' ^ , 

(4-  S  6  7  8 ) 

(b)  impregnated  nickel  sinter  type^  '  and  (c)  electro- 

(9  10  11  12 ) 

chemical  impregnation  ’  ’  .  The  last  method,  however,  is 

simple  in  operation  and  can,  in  theory,  achieve  high  loading  of 
active  material  in  a  single  impregnation  cycle  at  the  lowest  cost. 

In  the  electrochemical  impregnation  process,  electro- 
chemically  active  material  is  introduced  into  the  porous  nickel 
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plaqe  by  electrochemical  precipitation.  The  plaque,  which  serves 
as  a  cathode,  is  submerged  into  an  electrolysis  cell  containing 
an  aqueous  solution  of  cadmium  nitrate  at  a  temperature  near 
the  boiling  point  of  the  solution.  The  initial  pH  value  of  the 
Impregnation  solution  is  adjusted  to  around  3  by  titration  with 
nitric  acid.  (As  the  pH  value  of  the  bath  increases,  i.e., 
becomes  more  basic,  cadmium  hydroxide  begins  to  precipitate  in 
regions  other  than  within  the  pores  of  the  plaque,  and  active 
material  is  lost  in  the  bath,  which  is  not  what  v/e  want.) 

The  electrochemical  reaction  at  the  cathode  during  impre¬ 
gnation  could  either  be  the  liberation  of  hydrogen, 

2H2O  +  2e"  - ►  H2  +  20H"  ,  (1) 

or  the  reduction  of  nitrate  ion. 


(2) 


Both  reactions  produce  hydroxy  ions  which  then  either  diffuse 
away  from  the  electrode  porous  surface  or  react  with  the  cadmium 
ions  present  in  the  pores  to  form  cadmium  hydroxide.  The  cadmium 


H2O  +  NO^  +  2e  - NO2"  +20H 

•  t  •  •  • 

2H2O  +  NO^"  +  2e~  -  HNO2  +3OH 

t  •  •  •  • 

6H2O  +  :i0^"  +  8e" - NH_  +9OH 
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hydroxide  then  deposits  within  the  pores  of  the  porous  plaque. 

The  precipitation  reaction  iss 

Cd"^^  +  20H”  - »-  Cd(0H)2.  (3) 

There  is  the  possibility  of  forming  cadmium  according  to  the 
reduction  reactions 

Cd(0H)2  +  2e“  - *-  Cd  +  20H".  (4) 

This  reaction  is  presumably  prevented  from  taking  place  by 
hydroxy  ion  generated  on  the  electrode  surface  by  reaction  (1) 
or  (2),  so  that  there  is  no  chance  for  cadmium  ions  to  diffuse 
through  the  solution  onto  the  surface  and  be  reduced  to  cadmium 
metal . 


Following  impregnation,  the  electrode  is  cathodically 
polarized  in  a  25  potassium  hydroxide  solution  with  the  impre¬ 
gnated  plate  connected  as  the  cathode  and  subjected  to  low  direct 

2 

anodic  current  (0.25  to  0.75  arap/in  ).  The  solution  is  maintained 
at  about  100°C  for  a  period  of  from  15  to  45  minutes  to  convert 
cadmium  hydroxide  to  cadmium  in  the  pores  of  the  plaque  by 
reaction  (4)  or  to  remove  any  nitrate  ion  or  other  undesirable 
reduction  products.  The  active  material  of  the  electrode  is  pre¬ 
sumed  to  be  both  cadmium  hydroxide  and  cadmium. 
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In  this  work,  the  electrochemical  impregnation  process 
for  the  fabrication  of  the  cadmium  electrode  was  examined  both 
experimentally  and  theoretically. 

Experimentally,  a  set  of  current-density- vs.-time  transients 
at  constant  applied  potential  were  obtained  for  a  nickel  micro¬ 
electrode.  The  current  density  v/as  due  to  the  occurrence  of  the 
electrode  reaction  (2)  and  precipitation  reaction  (3)  presented 
earlier.  For  the  purpose  of  identifying  the  electrochemical 
reactions,  another  set  of  current-density -vs.-time  transients  were 
obtained  using  an  electrolyte  which  contained  no  cadmium  ion. 

In  this  case,  the  precipitation  reaction  (3)  no  longer  took 
place.  The  reactions  at  the  electrode  were  the  only  reactions 
taking  place. 

The  transient  current  responsed  in  an  electrolyte  which 
contained  cadmium  ions  differed  significantly  from  that  in  an 
electrolyte  which  did  not  contain  the  cadmium  ions.  This  diff¬ 
erence  is  presumed  to  be  due  to  the  precipitation  reaction. 

The  current  obtained  in  the  presence  of  cadmium  ions  in  the 
electrolyte  at  any  time  was  higher  than  that  obtained  at  the 
same  applied  potential  in  the  absence  of  cadmium  ions.  This 
result  is  explainable  by  the  fact  that  the  cadmium  ions  in  the 
solution  consumed  the  hydroxy  ions  which  were  produced  by  the 
heterogeneous  reaction  (2)  and,  as  a  result,  promoted  the  hetero¬ 
geneous  electrode  reactions. 
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The  experimental  current-density- vs,-time  curves  obtained 
in  the  potential-step  experiment  were  used  to  identify  both 
the  heterogeneous  electrochemical  and  the  homogeneous  preci¬ 
pitation  reaction  rate  expressions.  The  heterogeneous  rate  was 
identified  first  from  the  experimental  results  in  which  there 
was  no  cadmium  ion  in  the  electrolyte.  The  precipitation  rate 
was  then  identified  from  the  experimental  results  in  which  the 
cadmium  ions  were  present  in  the  electrolyte. 

A  set  of  differential  equations  which  described  the  trans¬ 
port  process  in  the  cadmium-ion-free  solution  was  solved  in  terms 
of  two  unknown  constants.  The  unknown  constants,  related  to  the 
heterogeneous  electrode  reactions,  were  identified  by  matching 
the  theoretical  results  with  the  results  of  the  experiment.  The 
precipitation  reaction  rate  constant  was  then  obtained  by  a 
similar  procedure. 
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CHAPTER  2 


LITERATURE  REVIEW 

Impregnation  of  nickel  sinters  with  cadmium  hydroxide, 
followed  by  electrochemical  reduction  to  cadmium,  is  the  most 
economical  method  by  which  to  fabricate  cadmium  electrodes. 

Almost  all  of  the  research  works  concerning  this  subject  may 
be  divided  into  three  types:  (a)  manufacturing  techniques  de¬ 
velopment,  (b)  impregnation  process  study  arid  (c)  charge-discharge 
study . 


During  recent  years,  mainufacturing  techniques,  especially 
those  involving  electrochemical  impregnation,  have  employed 
widely  varying  operating  conditions^ There  are 
several  parameters  that  affect  the  electrochemical  impregnation 
process  such  as: 

1)  Temperature:  Beauchamp  emphasized  the  fact  that  high 

impregnation  bath  temperature  keeps  the  size  of  cadmium 
hydroxide  crystals  small,  and  thereby  obtaining  a  more  active 
electrode . 

2)  pH  value:  In  Beauchamp's  process  ,  the  pH  value  was 

closely  controlled  by  buffering  the  cadmium  nitrate  solution 
with  sodium  nitrite.  Pickett  '  in  a  similar  process, 
maintained  the  pH  value  of  the  cadmium  nitrate  bath  between 
3-5  by  titrating  the  solution  with  proper  amounx  of  nitric 
acid . 
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3)  Current  density:  The  current  density,  based  on  apparent 

(14) 

electrode  surface,  is  another  important  variable.  Beauchamp 
used  D.  C.  current  of  4-8  amps/in^  while  Bulan  used 

2 

alternating  current  of  1.2-1. 6  amps/in  and  claims  that  the 
electrode  impregnated  by  alterning  current  technique  retained 
more  of  its  capacity  after  cycling  as  compared  to  electrodes 
fabricated  by  other  methods.  Pickett^ presents  another 
technique  similar  to  the  above  and  claims  that  his  technique 
yields  higher  loading  cadmium  electrodes  and  achieves  about 
85^  utilization  of  the  active  material. 

The  electrochemistry  of  the  Cd/Cd(0H)2  electrode  in  con¬ 
centrated  alkaline  solution,  typically  KOH  as  in  the  Ni-Cd 
battery,  has  been  studied  extensively  both  theoretically  and 
experimentally  ^  ^ .  Bennion,  et  al.^^'^^  in  their 

theoretical  development,  present  two  models,  the  solution- 
diffusion  model  and  film  model,  in  considering  the  failure 
mechanism  in  the  performance  of  secondary  batteries  using  metal/ 
metal  salt  couples.  They  conclude  that  one  mode  of  failure  of 
the  metal/metal  oxide  porous  electrode  is  the  blooking  of  the 

pore  or  the  complete  covering  of  surface  by  product  crystallities . 

(18) 

Falk'  obtained  X-ray  diffraction  patterns  from  electrodes 
submerged  in  electrolyte  during  charge  and  discharge  by  means 
of  a  special  test  cell.  Ke  claims  that,  during  charge,  Cd(CH)2 

is  successfully  transformed  into  Cd  metal  in  the  cadmium  electrode, 
but  even  after  a  strong  overcharge,  this  transformation  is  not 
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complete,  as  indicated  by  the  X-ray  pattern  at  this  stage. 

Furthermore,  Cd(OH),  is  the  end  product  during  discharge  and 

there  is  no  evidence  of  the  presence  of  a  CdO  film  which  is 

believed  to  be  responsible  for  the  passivation  of  the  electrode. 

( 19 ) 

Gross  and  Clocking  ^  summarized  the  results  of  recent  research 
by  characterizing  the  negative  cadmium  electrode  in  a  nickel- 
cadmium  cell  and  evaluating  some  of  the  published  results  about 
which  researchers  still  disagree  among  themselves,  such  as  the 
mechanism  for  passivation,  which  is  an  important  phenomena  that 
affects  the  capacity  of  the  nickel-cadmium  cell. 

In  comparison  to  the  extensive  studies  on  the  electro¬ 
chemistry  of  the  Cd/Cd(0H)2  electrode  in  KOH  solution,  there 
is  little  information  available  concerning  the  electrochemical 
reactions  taking  place  during  the  electrochemical  impregnation 
process  which  uses  cadmium  nitrate  solution  as  electrolyte. 

The  only  research  works  known  to  the  author  are  the  electro- 

(21  22 ) 

cnemical  impregnation  studies  using  cycling  voltammemry  ’  . 

2  2 
They  used  a  0.013  cm  nickel  microelectrode  and  a  0.064  cm 

cadmium  microelectrode.  Some  important  conclusions  of  their 

studies  are  summarized  below: 

(1)  There  was  a  large  reduction  wave  at  -0.d3V  vs  SGE( Saturated 
Calomel  Electrode)  which  appeared  only  during  the  initial 
negative  scan.  This  is  believed  to  be  the  reduction  of  a 

solution  species  to  liberate  0H~  ion  in  the  presence  cf 

+2 

Cd  to  precipitate  cadmium  hydroxide  on  the  electrode 


surface . 
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(2)  A  second  reduction  process  was  observed  at  -0.62V  vs  SGE 
during  the  reverse  positive  scan.  This  is  due  to  the  reductive 
formation  of  metallic  cadmium  from  cadmium  oxide  or  hydroxide 
that  was  deposited  on  the  electrode  surface  during  the  pre¬ 
vious  forward  scan. 

(3)  Some  processes  (or  one  of  the  processes)  during  the  first 
cycle  effectively  passivated  the  electrode  to  further  electro¬ 
activity.  The  formation  of  the  gray  material  at  -O.63V  vs 
SCE,  discussed  in  (1),  was  the  most  probable  reason  for  the 
passivation  of  the  electrode. 
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CHAPTER  3 

THEORETICAL  APPROACH 

3  •  1  General 

In  voltammetry  at  constant  voltage,  the  current  through 
the  electrolytic  cell  is  measured  as  a  function  of  time  while 
the  potential  of  the  polarized  electrode  is  held  at  a  constant 
value.  The  current  is  affected  by  both  the  transfer  process  and 
the  electrochemical  and/or  chemical  reactions. 

A  plane  electrode  is  immersed  in  a  solution  which  contains 
a  substance  0.  At  the  imposed  potential,  0  is  reduced  to  another 
substance  R  at  the  electrode  surface,  according  to  the  following 
electrochemical  reactions 

0  +  ne~  - ►  R  (5) 

The  current- density- vs.-time  curve  obtained  at  a  constant 
voltage  depends  on  both  the  kinetics  of  the  electrochemical 
reaction  and  the  rate  of  mass  transfer. 

Under  our  experimental  conditions,  diffusion  is  the  only 
mode  of  mass  transfer.  The  other  two  modes  of  mass  transfer, 
namely,  migration  and  convection,  are  insignificant.  Migration 

I 
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is  eliminated  by  the  presence  of  a  large  excess  of  supporting 
electrolyte  which  reduces  the  potential  field.  Convection  is 
avoided  due  to  the  conditions  that  (1)  the  solution  is  not 
stirred  and  (2)  the  duration  of  electrolysis  is  short,  so  that 
the  density  change  has  not  yet  become  a  significant  factor  and 
natural  convection  has  not  taken  place. 


The  dimension  of  the  cell,  as  compared  to  the  microelectrode 
used,  is  so  large  that  the  solution  may  be  regarded  as  extending 
to  infinity;  i.e.,  the  diffusion  process  is  semi-infinite. 
According  to  Pick's  first  law,  the  flux  of  the  substance  0  at 
a  distance  x  from  the  electrode  and  in  the  direction  perpen¬ 
dicular  to  the  electrode  is: 

'^CAx,t) 

'•o  '  -  °o-t - 

Ox  • 

Conservation  of  the  species  0  leads  to  the  following  differential 
equation ; 

'aCn(x,t)  ^ 

t  °  ,  (?) 


where  Cq  is  the  concentration  of  substance  0,  which  is  a  function 
of  x  and  t.  Dg  is  the  diffusion  coefficient  of  substance  0, 
and  is  assumed  to  be  a  constant. 
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When  the  substance  0  is  reduced  at  a  very  high  rate,  by 
imposing  a  sufficiently  large  potential  step,  such  that  the 
concentration  of  this  substance  at  the  electrode  surface  is 
reduced  to  zero  soon  after  the  start  of  electrolysis ,  the 
boundary  condition  at  the  electrode/electrolyte  interface  for 
equation  (?)  is  CQ(0,t)=0  for  t>0.  Initially  before  the  electro¬ 
lysis,  the  solution  is  homogeneous,  and  the  concentration  Cq 
at  t=0  is  uniform  with  the  bulk  value.  The  other  boundary 
condition  is  that  Cq  approaches  its  bulk  value  when  x  is 
sufficiently  large. 

The  electrolysis  current  is  equal  to  the  product  of  the 
charge  involved  in  the  reduction  of  one  mole  of  substance  0  by 
the  flux  of  this  substance  at  the  electrode  surface.  Thus  the 
current  density  (current  per  unit  electrode  surface  area)  is: 

1=  nFNQ(0,t),  (8) 

where  n  is  the  number  of  electrons  involved  in  the  reduction 
reaction  .of  substance  0,  F  is  the  Faraday  constant,  and  NQ(0,t) 
is  the  flux  of  substance  0  for  x=0. 

The  solution  of  the  above  problem  is^^^^  : 

1=  nFD  i  cl  - -  (9) 

t 

where  Cq  is  the  bulk  value  of  substance  0. 
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The  current  density  represented  by  equation  (9)  is  called 
the  limiting  current  density  of  electrochemical  reaction  (5)^ 

This  is  the  maximum  current  density  which  can  be  achieved  by 
this  reaction.  It  was  found  that  if  nitrate  ion  is  the  only 
reducible  substance,  the  current  density  predicted  by  equation 
(9)  is  considerably  lower  than  the  experimental  value.  This 
contradition  is  resolved  if  one  realizes  that  under  the  experi¬ 
mental  conditions  water  is  also  reduced. 

At  any  other  less-negative  potential  the  current  will 
depend  on  the  electrode  kinetics.  Now  the  condition  that  the 
reducible  substance  concentration  becomes  zero  after  the  appli¬ 
cation  of  the  potential  step  has  to  be  replaced  by  a  kinetics 
equation,  which  is  a  function  of  potential  and  the  concentrations 
of  surface  reactant  and  product.  No  general  solution  is  available 

except  for  the  case  when  the  reaction  rate  is  linearly  dependent 

( 21 ) 

on  tne  concentrations  of  the  surface  reactant  and  product  ^  . 

3.2  Description  of  the  Model 

The  real  porous  electrode  is  formed  by  sintering  nickel 
powder  on  a  nickel  grid  to  form  a  plaque.  The  plaque  serves  as 
the  cathode,  which  is  positioned  in  a  cadmium  nitrate  solution 
with  a  pH  value  between  3  to  5*  The  passage  of  the  current 
produces  hydroxy  ions  which  cause  the  cadmium  to  be  deposited 
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into  the  porous  pores.  The  transport  phenomena  associated 
with  electrochemical  reactions  in  a  porous  medium  is  very 
complicated.  The  approach  here  is  to  study  the  intrinsic 
reactions  which  take  place  inside  the  pores  in  a  simple 
flat  electrode.  The  plate  electrode  is  submerged  in  a 
semi-infinite  pool  of  electrolyte  so  that  the  mass  transfer 
problem  can  be  treated  as  a  one-dimensional  problem  in 
the  x-direction  (perpendicular  to  the  electrode  surface). 

In  the  solution  side,  there  is  a  thin  double-charge  layer 
near  the  electrode.  The  thickness  of  the  double  layer  is 

O 

on  the  order  of  10  to  100  A.  It  is  too  thin  to  be  probed 
adequately,  and  the  theory  of  the  diffuse  layer  is  a 
microscopic,  rather  than  a  macroscopic,  phenomena.  The 
mass  transfer  region  to  be  considered  here  is  thus  located 
outside  of  this  double  layer. 

The  electrodeposition  process  operates  on  the  principle 
that  the  hydroxy  ion,  which  is  generaxed  by  the  reduction 
process,  coprecipitates  with  the  cadmium  ion  in  the  solution 
to  form  cadmium  hydroxide  crystalloid.  The  approach  taken 
here  to  unravel  the  complicated  sequence  of  reactions  is 
to  form  a  model  which  describes  the  transport  processes 
and  reactions.  The  model  will  first  be  used  to  analyze 
the  experimental  results  obtained  under  the  condition  of 
no  pricipitation  reaction,  by  using  an  electrolytic  solution 
which  is  free  of  cadmium  ions.  The  heterogeneous  reduction 
reaction  will  be  determined.  The  transport  and  reacxion 
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problems  will  then  be  examined  in  the  case  v/hen  the  solution 
contains  cadmium  ions.  The  rate  of  the  precipitation 
reaction  will  then  be  determined  by  matching  the  model's 
prediction  and  the  experimental  results. 


3 . 3  Formulation  of  Transport  Process  Equations 


Three  ionic  species  to  be  considered  in  the  electro- 

-  -  +2  • 

chemical  system  are  NO^,  On  and  Cd  ions.  In  the  absence 
of  migration  and  convection,  the  flux  of  each  species  is 
expressed 


h’i  -  -  D, 


,  i=l,^.3  , 


where 


is  t;;e  flux  of  species  i  (moles/cm^-sec ) , 

is  the  diffusion  coefficient  of  species  i  (cm  /sec), 

is  the  concentration  of  species  i  (raoles/cm^), 

which  is  a  function  of  both  x  and  t,  and  i=l,2,3 

-  -  +2 

corresponding  to  NO^,  CH  and  Cd  ions. 

(24) 

The  equation  of  continuity  of  each  species  is 


I2x  = 

"St 


d 


i=l,2,3  , 


v/here  is  the  rate  of  production  of  species  i 
( moles/cm^-sec )  by  precipitation  reaction,  which  occurs 
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only  in  the  presence  of  cadrr.ium  ions. 


Substitution  of  equation  (10)  into  (11)  gives 


2£i 

?>t 


2 

^  C. 


=  Di' 


2):. 


2 


i=l,2,3 


(12) 


It  is  assumed  that  the  precipitation  reaction  (3) 
can  be  treated  as  a  homogeneous  reaction  and  is  linearly 

proportional  to  the  degree  of  supersaturation.  Therefore, 

the  consumption  rate  of  cadmium  ion  due  to  this  reaction 

is  ; 


R 


3 


(13) 


where  k  is  the  homogeneous  reaction  rate  constant  (sec"^), 
K sp  is  the  solubility  constant  of  Cd(0H)2 

and  has  the  value  of  2x10”^^  (moles^/cm^) 

3y  the  stoichiometric  relation  of  equation  (3)1  '^^e  con¬ 
sumption  rate  of  hydroxy  ion  due  to  the  same  reaction  is  : 


rt2  -  2  R3  - 


-2k 


c.  -■ 


(C^)^ 


(14) 


Substitution  of  equations  (I3)  and  (14)  into  .equation  (12' 
leads  to  the  following  two  differential  equations  i 


Dt 


5^-  2  k 


Ksp 


N 


(C^)' 


V  15) 


2  (' 

k  C- 

I  3 


/  p  '2 


For  N'O"  ion, 


=  0. 


Equation  (12)  for  NO^  ion  may  be  reduced  to 


■2)^1 

- —  =  D  - ^ 

•at  1x2 

O  X 


Only  two  of  the  three  differential  equations  need 

to  be  solved  for  the  concentration  profiles,  due  to  the 

(24) 

electroneutrality  condition  : 


203-0^-02=0 


Substitution  of  equation  (18)  into  equation  (I5) 


gives  s 


^2 

at 


Equations  (1?)  and  (19)  are  solved  simultaneously 

to  obtain  the  concentration  profiles  of  0^  and  O2,  i,e., 

-  .  +2 

..■03  and  OH  ions.  The  0^  profile  can  then  be  obtained  oy 
equation  (18)  from  the  knov/n  profile^  of  and  C2. 


A  set  of  differential  equations,  wnich  are  simpler 


18 


than  equations  (17)  and  (19)i  arsi  obtained  when  the  elec 
trolyte  contains  no  cadmium  ion.  In  the  absence  of  cad¬ 
mium  ion,  the  precipitation  reaction  (3)  can  not  occur; 
thus,  terms  are  all  equal  to  zero.  Equation  (12)  is 
then  reduced  to: 


•St 


i=l,2  , 


(20) 


in  which  i=l,2  corresponds 


to  NO^ 


and  OH 


ions . 


Equation  (20)  is  written  for  and  C^: 


■s^i 

■St 


■St 


(21) 


(22) 


The  other  condition,  the  concentrations  of  various 
ions  at  the  electrode/electrolyte  interface  after  the 
start  of  the  current  flow,  has  to  be  described  by  the 
electrode  kinetics  of  the  reduction  reaction  of  N0“  ions 


3.^  Electrode  Kinetics 


The  possible  reduction  reactions  are  numerous,  as 
shown  in  equation  (2).  Since  in  e.-ch  reaction,  OH"  ion 
is  generated,  we  will  represent  the  production  of  0H“ 
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ions  as  the  following  electrochemical  reacrior. : 

+  NO^  T  2e“ - ►  +  30H“.  (23) 

Because  there  are  two  electrons  involved  in  the  reaction, 
it  is  likely  not  an  elementary  step.  Let's  assume  the 
reaction  comprises  two  elementary  steps,  each  step  involving 
the  transfer  of  one  electron,  and  a  dissociation  reaction: 


(slow) ,  (2La) 

(fast),  (2Lb) 

(fast),  (25) 


where  and  cathodic  reaction  constants,  and 

and  K^2  anodic  reaction  constants,  respectively.  is 
the  dissociation  reaction  constant.  unstable 

intermediate  ion  which  immediately  enters  the  next  reaction  (2Lb) 
The  rates  of  these  elementary  steps  should  always  be  propor¬ 
tional  to  each  other.  For  the  elementary  steps  listed  above, 
reaction  (2Lb)  should  occur  once  every  time  reaction  (2da)  occurs 


Let  I.  and  I,  denote  the  cathodic  current  densities 

X  K. 

of  reactions  (2La)  and  (2Lb),  and  r^  and  r^  denote  che 
reaction  rates  of  reactions  (24a)  and  (24b),  respectively. 
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Then,  "ne  relationship  between  ,  I.  and  surface  ion 


concentrations  of  I'iCi,  (CH)^  ,  v/hich  app 

( 2  ) 

step  (24a),  may  be  v/ritten  as  ; 


ear  in  eiementar; 


.  -  i. 

'l"  nF 


(n=i: 


=  -41  [;=Nc;]|??o:-:)-]  f- 


[1-^1  :f 


..  r.o  1  r„o  *1  r  £ 

cl  r^OZ  rl^o  ^ 

L  JJ  w  ''  L  *'• 


r 


Similarly,  for  reaction  (24b),  one  has: 


(n=l} 


Kc2  [=(0K)-] 


where  and  are  symmetry  factors  representing  the  fractions 
of  applied  potential  V,  which  promotes  the  cathodic  reaction. 

an^l  ^ru-  sne  the  surface  ion  concentrations  of 

^  U** 

:<0y  (OH)^  and  OK"  on  the  electrode  surface. 


Since  reactions  (24a)  and  (24b)  occur  at  the  same  rate, 
we  have : 


1=  =  2I2  =  21.  , 


2 
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where  I  is  overall  curren'  density.  It  is  further  assur.ed 
that  reaction  (24b)  is  fast,  and  is  essentially  in  equilibrium. 
From  equation  (27), 


^c2 


and  = 
r 


=r  =  overall  reaction  rate. 


Substituting  equation  (28)  into  equation  (26)  and 
combining  equation  (29)  gives; 


_  ^al^a2  ('  ,.o 


(2-/3^)? 


From  reaction  (25)  one  may  write 


*^d  = 


['-:%-]  pH^cj 


Substitution  of  the  above  expression  into  equation  (32)  leads 


22 


23 


ions  and  electrolyte  containing  no  cadmiuni  ions,  ~vo 
and  four  boundary  conditions  are  needed  to  solve  the 
equations  (1?)  and  (19)  or  (21)  and  (22). 


initial 
dif f er en 


At  t=0,  the  concentration  of  each  ion  is  equal  to  its 
bulk  value  at  each  point  along  the  x-direction.  Ihe  initial 
conditions  are  as  follows: 

1.  at  t=0,  C^=  for  all  x, 

2.  at  t=0,  €2=  C2  for  all  x, 

where  the  superscript  "b"  denotes  the  bulk  condition, 


The  first  boundary  condition  is  the  consequence  of 
the  fact  that  the  flux  at  the  interface  has  to  be  equal  to 
the  rate  of  heterogeneous  electrode  reaction,  i.e.: 

where  superscript  "0"  denotes  the  surface  condition.  To 
satisfy  the  stoichiometric  relation  of  equation  (23),  one 
has,  at  x=0 : 


X 

f  \ 

^x 

1 

0 

II 

^2  -v 

^  J 

This  is  the  second  boundary  condixion.  The  cx;.er  r.vc 
boundary  conditions  are; 
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3.  at  x=oo  ,  C^=  C°, 

•L. 

4.  at  x=oo  ,  C^=  C^. 

That  is,  the  concentrations  far  away  from  the  electrode/ 
electrolyte  interface  do  not  change. 


The  flux  in  equation  (32)  is  proportional  to  the  curren 
density  as  a  function  of  time: 
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CHAPTER  4 

EXPERIMENTAL  WORKS 

A  potential-step  was  generated  from  a  function 
generator  and  was  applied  to  the  working  electrode/elec¬ 
trolyte  interface  in  the  electrochemical  cell.  The  ca¬ 
thodic  current  was  thus  sampled  and  stored  in  memory  at 
a  constant  time  interval  by  a  microcomputer.  The  resulting 
Gurrent/time  transient  was  then  plotted  by  a  x-y  recor¬ 
der  or  printed  out  by  a  line  printer  by  recalling  these 
current  data  from  the  memory  of  the  microcomputer. 

4. 1  The  Electrochemical  Cell 


The  current  transient  under  a  potential  step  con¬ 
dition  was  obtained  in  a  three-electrode  cell.  The  cell 
and  electrodes  configuration  is  sho-wn  in  Fig-ure  1.  The 
working  electrode  was  a  nickel  microelectrode.  It  was  made 

of  a  nickel  wire  which  was  pressure-fitted  into  a  teflon 

2 

cylinder.  The  flat,  exposed  surface  was  0.013  cm  .  A 
saturated  calomel  electrode  (SCE)  and  a  cadmium  bar  were 
used  as  the  reference  and  the  counter  electrodes,  res¬ 
pectively  . 


The  cell  was  filled  with  electrolyte  of  either 
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cadmium  nitrate  solution  or  potassium  nitrate  solution. 

The  pK  value  of  the  solution  v/as  controlled.  Potassium 
chloride  was  used  as  the  supporting  elecxrolyte.  The  electro¬ 
lytic  solution  was  not  stirred  during  the  experiment,  so  that 
the  conditions  of  serai-infinite  linear  diffusion  were  maintained. 

Two  experimental  parameters  were  varied: 

1)  the  ion  concentration  of  the  uolurions  and 

2)  the  magnitude  of  the  potential-step  applied  zo  the 
interface. 

4,2  The  Instruments  and  Electrical  Setup 

The- block  diagram  for  the  experiment  setup  is  sho'wn  in 
Figure  2.  The  instruments  sho'wn  in  Figure  2  contain  the 
following  equipment: 

a)  A  Princeton  Applied  Research  (?AF)  Model  Potentiostat/ 
C-alvanostat :  In  the  experiment,  the  operating  mode  was  set 
at  CONTROL  E.  ,  which  means  that  the  instrument  f-unctioned 
as  a  potentiostat.  The  current  was  measured  while  keeping 
the  potential  of  the  working  electrode  constant.  This  was 
accomplished  by  setting  the  counter  electrode  to  the 
required  level  so  that  tne  v/orking  electrode  potential  was 
at  a  programmed  potential  ■■/itn  respect  to  the  reference 
electrode.  The  instrument  is  prcviied  v/ich  a  model  l"o 
Current- to-7oltage  Converter  which  is  the  "basic"  plug-in 
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module.  The  M  176  provides  a  do  vol-age  output  -./hi oh  is 
proportional  to  the  current.  The  cell  current  can  be  eitner 
displayed-  on  the  panel,  by  a  x-y  recorder,  or  monitored  by 
the  microcomputer.  Connection  to  the  cell  was  made  through 
the  external  cable »  The  counter  electrode  v/as  ccrn.ected  to 
the  red  clip,  the  working  electrode  to  the  green  clip  and 
the  reference  electrode  to  the  Electrometer  Probe  which  has 
a  very  high  impedance,  thereby  insuring  us  that  the  current 
was  flowing  to  the  reference  electrode. 

b)  A  PAP  Model  175  Function  Generator:  This  is  a  programmable 
waveform  generator  which  has  two  operating  modes,  SV/EE?  and 
PULSE.  The  former  generates' a  sequence  of  triar.gialar  waves, 
while  the  latter  generates  a  sequence  of  step  function.  For 
the  potential  step  experiment,  the  operating  mode  was  set-  a 
PULSE  mode  and  the  pulse  width  selector  at  STEP  position. 
The  magnitude  of  the  step* function  was  set  by  the  setting  o 
the  3  potential  (upper  limit)  in  the  PCTSMTIAL  section  on 
the  front  panel  of  this  instrument. 

c)  An  Intel  CPU  8080A  Bases  Tv'licrocom'Duter  ( c-lK  PA.'<.):  The  micro 
computer  was  the  central  part  of  the  experimental  setup. 

The  functions  of  the  microcomputer  in  the  experiment  were: 

(1)  to  trigger  the  Function"' Cenerator  which  activated  the 
potential  step  presented  earlier. 

(2)  -to  sample  the  current  outputs  from  tne  M  17t  Converter 

-  •  -  and  store--them  irv- the- memory .  — 

(3)  to  convert  these  digital  data  back  into  analog  signals 
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and  then  plot  them  on  the  x-y  recorder. 

(ir)  to  print  the  digital  data  on  the  line  printer. 

(5)  to  store  these  digital  data  in  the  tapes. 

The  peripheral  equipment  of  the  microcomputer  include  a 
CRT  (TV  screen),  the  keyboard,  a  cassette  recorder  and  a 
line  printer. 

d)  A  Houston  Model  RE  007^  x-y  Recorder;  The  Recorder  receives 
the  analog  outputs  from  the  DAC  (Digital-Analog  Converter) 
and  plots  them  on  the  graph  paper. 

4.3  The  Computer  Program 

A  computer  program  'vvritten  in  BASIC  language  was  used  xo 
carry  out  the  experiment.  The  program  consists  of  a  main  program 
and  a  machine-language  subroutine  called  "MUG".  This  subroutine 
includes  an  execution  statement  which  generates  a  trigger  signal 
to  activate  the  Function  Generator.  A  potential  step  was  xhen 
immediately  applied  to  the  working  electrode,  and  thereafter 
the  computer  started  to  sample  and  store  the  current  output  at 
a  fixed  time  interval.  The  analog  signal  was  converted  into  a 
digital  datum  by  a  12  bits  ADC  (Analog-Digital  Converter)  and 
then  stored  in  the  memory.  When  the  number  of  samples  reached 
a  preset  value,  the  sampling  routine  was  terminated.  The  digital 
signals  stored  in  the  memory  were  binary  numbers.  These  data 
v/ere  converted  into  decimal  numbers  and  then  converted  into 
analog  signals,,  and  finally  plotted  on  the  x-y  recorder  and 
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printed  on  the  line  printer.  This  was  executed  by  the  main 
program.  The  main  program  and  the  "MUG"  subroutine  are 
presented  fully  in  Appendix  A1  and  A2 . 

4.4  Experimental  Conditions  and  Procedures 


In  the  experiment,  the  initial  pH  value  of  the  eieccro- 
lyte  was  maintained  at  3*0,  which  was  mhe  condition  used  in  the 
actual  electrochemical  impregnation  process.  This  v/as 
accomplished  by  titrating  the  electrolytic  solution  with  a 
concentrated  HCl  solution.  The  concentrations  of  nitrate 
ion  used  were  0.005M,  0.0025rvi  and  0.00125M  prepared  from  either 
the  or  Cd(NO^)  solution.  The  cadmium  ion  concentrations 

were  0.0025M,  0,00125M  and  O.OOQ625M,  All  the  electrolytic 
solutions  contained  0,5M  KCl  as  the  supporting  electrolyte. 

It  was  .found  later  that  the  electrolysis  of  water  made 
significant  contribution  to  the  total  current,  because  the 
current  obtained  was  significantly  higher  than  the  limiting 
current.  In  order  to  obtain  the  current  due  only  to  the  reduction 
of  nitrate  ions,  several  ..r-uns  with  solutions  which  contained  no  . 
nitrate  ion  were  conducted.  The  solutions  used  for  this  purpose 
v/ere  either  water  or  C.GC23M  CdGl^  solution.  Bctn  solutions 
contained  C.yM  KCl  as  tne  supporting  electrolyte.  This  current 
was  subtracted  from  the  current  obtained  in  tne  oresence  of 
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nitrate  ions  in  the  solution.  The  resulting  current  is  the 
faradic  current  of  the  reduction  of  nitrate  ions. 

The  magrxitudes  of  the  potential-step  used  in  this 
study  were  -'0V40V,  -0 .60^ 'and  ' -0.80V  versus  the  equilibrium 
potential  which  was  at  approxix-nately  -0.35^'^  vs  SCE, 

The  experimental  procedures  are  listed  blow; 

(1)  The  electrical  circuit  was  set  up  as  shown  in  Figure  2, 
except  that  the  cell  was  disconnected  from  the  Potentiostat/ 
Galvanostat  by  setting  the  cell  selector  on  the  i.nstrument 
to  the  OFF  position.  A  proper  current  range  was  then 
chosen  from  the  Current  Range  Switch  to  make  sure  that 

the  current  output  would  not  be  overloaded.  The  Function 
Generator  was  initialized  by  depressing  the  INITIAL 
control  pushbuttom.  The  initial  potential  was  set  at 
zero  volts. 

(2)  A  known  volume  of  electrolytic  solution  and  an  equal 
volume  of  supporting  electrolyte,  KCl,  were  added  to  the 
cell  and  titrated  with  HCl  solution  to  a  pH  value  of  3*0. 

The  solution  was  then  deaerated  by  bubbling  purified 
nitrogen  through  the  stirred  solution  for  about  ten 
minutes.  The  gas  continued  to  pass  above  the  solution 
during  tne  experiment. 

(3)  After  the  working  electrode  surface  was  polisned  by  a 

3/0  alumina  paper,  it  was  rinsed  thoroughly  with  distilled 
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and  deionized  water,  and  then  was  positioned  in  the  cell. 

The  electrodes  (the  working  electrode,  counter  electrode 
and  reference  electrode)  were  then  connected  to  the 
Potentiostat/Galvanostat  as  described  in  section/;. 3* 

After  the  electrodes  were  correctly  connected,  the  cell 
selector  was  switched  to  the  EXTERNAL  CELL  position.  The 
rest  potential  of  the  solution  could  be  determined  by 

turning  the  knob  in  the  APPLIED  POTENT I AL/CURRENT  section 
of  the  front-panel  of  the  Potentiosuat/Galvanostat  so 
that  the  current  reading  displyed  on  the  front  panel 
meter  was  zero.  The  applied  potential  was  then  the  rest 
potential. 

(4)  The  upper  limit  (B  potential)  on  the  Function  Generator 

was  set.  The  value  of  the  -3  potential  plus  the  rest 
potential  was  the  total  potential  applied  to  the  electror 
chemical  cell.  . 

(5)  Lhe  experiment  was  initiated  by  running  the  computer 
program.  The  current/time  data  were  stored  in  the  computer’s 
memory  and  then  the  results  were  converted  to  line 

printer  output. 

4. 5  Experimental  Results 

Two  sets,  a  total  of  twenty-two  runs,  were  made.  The 
first  set  of'  runs  was  conducted  in  solutions  which  contained 
no  cadmium  ion.  The  second  set  of  runs  was  made  in  solutions 
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which  contained  cadmium  ions.  For  each  set  of  runs,  clanK  runs 
were  made  in  the  solutions  without  the  presence  of  nitrate 
ion  for  the  purpose  of  determining  the  background  current 
due'to  the  water  decomposition  reaction.  The  solution  for 
blank  runs  for-the  first  set  was. water,  while  the  solution 
for  blank  runs  for  the  second  set  was  made  by  using  the 
O.OO25M  GdCl2  solution. 


The  current-density-vs .-time  data  printed  by  the  H  14 
Line  Printer  are  sho’wn  in  Table  3.1  through  Table  B.22  in 
appendix  3.  The  measured  current  was  converted  to  the  current 
density  which  appears  in  those  tables.  This  is  calculated 
by  dividing  the  currant  by  the  electrode  surface  area 
0.013cm^. 


Tables  3.1  through  3.9  are  the  first  data  set  which  used 
potassium  nitrate  as  the  electrolytic  solution,  while  Tables 
3.10  through  3,12  are  the  background  currents  of  water  used 
to  subtract  from  the  first  data  set. 


Tables  3.13  through  3.21  are  the  seccr.d  data  set  which 
used  cadmium  nitrate  as  the  electrolytic  solution,  while 
Table  3.22' is  the  background  current  of  water  in  t.ce  presence 
of  cad.mium  ions. 


A  series  of  current-density-vs . -ti.me  curves  witn  various 
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electrolytic  solution  concentrations  and  applied  overpotentials 
are  shown  in  Figures  3  through  6.  In  these  figures,  the 
current  due  to  the  electrolysis  of  water  has  been  eliminited 
by  using  the  routine  tabulated  in  Table  1. 

Figures  3  through  5  are  the  current-densivy-vs . -time 
curves  obtained  in  the  solutions  which  did  not  contain  cadmium 
ion.  Figure  3  shows  the  current-density-vs . -time  curves 
obtained  in  a  solution  containing  O.OO5M  nitrate  ion  and  at 
potentials  -0.40V,  -O.6OV  and  -O.SOV  from  the  equilibrium 
potential.  Figure  4  shows  the  current-der.sity-vs , -time  curves 
obtained  in  a  solution  containing  O.OO25M  nitrate  ion,  and  at 
potentials  -0.40V,  -O.6OV  and  -0.807  from  the  equilibrium 
potential.  Figure  5  shows  the  current-density-vs . -time  curves 
obtained  in  a  solution  containing  0.00125M  nitrate  ion  at  -0.40V, 
-O.6OV  and  -O.SOVfrom  the  equilibrium  potential.  Due  to  the 
lack  of  blank  runs  with  O.OOI25M  and  0.000625M  cadmium  ion  con¬ 
centrations  and  -O.6OV  and  -O.SOV  from  the  equilibrium  potential, 
only  one  current-density-vs . -time  curve  was  obtained  witn  0.005M 
nitrate  ion  and  0.0025M  cadmium  ion  concentraxion  at  the  poten¬ 
tial  of  -0.40V  from  the  equilibrium  potential.  This  is  shown 
in  Figure  6. 

6  Discussion  of  Exnerim.enxal  Results 


Some  characteristics  can  be  seen  from  the  curves  snov/n 
in  Figures  3  through  6. 
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Table  1:  Evaluation  of  Current-Density-vs . -Time  Data 
Plotted  in  Figures  3  Through  6. 


Figure  is 
plotted 

by  subtracting  the 
corresponding  current 
density  at  each  time 

Doint  in 

from  that 

of  in 

(1) 

. 

Table  3.1 

Table  B.  10 

Figure  3 

(2) 

Table  B.2 

Table  By  11 

(3) 

Table  3.3 

Table  3. 12 

(1) 

Table  3.4 

Table  3.10 

Figure  4 

(2) 

Table  B.5 

Table  3.11 

(3) 

Table  3.6 

Table  B.12 

I 


(1) 


Table  3.7 


Table  5.10 


Firstly,  the  current  density  decays  with  ti.te  in  a  ] 

hyperbolic  fashion.  This  is  the  result  of  depletion  of  the  | 

NO^  ions  near  the  electrode  surface.  Theoretically,  there  j 

is  an  initial  sharp  rise  of  current  associated  with  double-layer  j 

(27) 

Charging  ,  due  to  the  application  of  a  potential-step.  This  was  J 

not  observed  in  the  experiment,  however,  because  the  time  for  | 

charging  is  usually  very  short  as  compared  to  the  electrolysis  1 

time. 

i'l 

Secondly,  the  sharp  current  drops  in  the  beginning  indicate  j 

that  very  sharp  concentration  gradients  are  established  for  the 

I 

nitrate  ions  and  the  hydroxy  ions  as  soon  as  the  surface 

reaction  (23)  takes  place.  The  concentration  gradients  are  ! 

i 

gradually  decreased  due  to  the  diffusion  layer  extending  in  ^ 

the  direction  of  decreasing  concentration  gradients.  The  i 

result  is  that  the  decay  of  the  current  density  is  not  as  fast  | 

as  in  the  beginning.  * 


The  curves  for  the  case  of  potential  at  -0,80V  from  the 
equilibrium  potential  and  solutions  containing  0.0Q25:<1  and 
O.OOI25M  nitrate  ion,  as  sho-wn  in  the  upper-most  curves  in  Figures  - 
and  5»  respectively,  are  somewnat  different  from  tne  others. 

These  two  curves  do  not  follov/  t.te  same  trends  as  tneir  lower 
potential  counterparts.  It  is  suggested  that  tne  reactant, 
nitrate  ion,  near  the  electrode  surface,  is  used  up  in  a  very 
short  time  and  another  reaction  .must  be  taking  place.  This 


behavior  is  the  mosz  pronounced  as  sho'A/n  in  the  upper-most 
curve  in  Figure  5*  A  maximum  current  is  observed  at  isime 
equal  to  about  0.5  msec.  This  curve  thus  strongly  supports 
our  explanation  that  some  other  reaction  is  talcing  place  at 
that  time.  On  the  other  hand,  the  absence  of  this  behavior 
in  the  curve  for  the  case  of  O.OOpM  nitrate  ion  shO'*vn  in  the 
upper-most  curve  in  Figure  3  suggests  that  the  nimrate  ion 
concentration  neacr  the  electrode  surface  is  not  zero. 

Two  parameters  were  varied  in  the  experiment,  i.e., 
the  concentration  of  the  electrolyte  solutions  and  the  value 
of  the  potential  from  the  equilibrium  potential  (rest  potential). 
The  changing  of  either  the  bulk  concentration  of  the  electrolyte 
solutions  or  the  overpotential  that  is  applied  to  the  working 
electrode  will  change  the  value  of  the  current  density.  For 
a  given  bulk  concentration,  the  increase  in  the  applied  over- 
potential  increases  the  electrochemical  reaction  rate  on  the 
electrode  surface.  The  effect  of  the  applied  overpotential  on 
the  current  density  can  also  be  observed  in  Figures  3  through  5> 
In  each  of  these  figures,  a  plot  of  the  limiting  current - 
density- vs .-time  curve  is  also  shown.  The.  limiting  current 
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zero  at  any  tine.  After  examining  Figures  3  through  5i  find 
thatj  for  the  case  of  -0.30Y,  some  of  'he  values  of  current 
density  exceed  the  limiting  currer.o  density.  This  phenomenon 
becomes  more  significanx  when  the  bulk  concentraxion  of  poxa- 
ssiura  nitrate  is  decreased.  The  reason  for  this  phenomenon 
probably  is  that  at  the  high  applied  overpotentiai ,  many 
reactions  can  subsequently  be  expected  to  take  place.  The 
maximum  current  appearing  in  xhe  curve  for  the  case  of  -O.cOV 
is  a  resulx  of  this  behavior. 

IVhen  xhe  applied  overpotentiai  is  fixed,  the  hetero¬ 
geneous  reaction  rate  constants  in  boxh  the  cathodic  and  anodic 
directions,  and  K2 ,  are  fixed.  The  supply  of  the  nixrate 
ions  on  the  electrode  surface  comes  only  in  the  way  cf  diffusion 
of  the  nitrate  ions  from  the  bulk  soluxion.  Increasing  the 
bulk  concentration  will  increase  the  raxe  of  nitraxe  ion  supply 
to  the  electrode  surface,  thusly  increasing  the  current  densixy. 
Figures  7  xhrough  9  show  the  effect  of  the  bulk  concenxraxion 
of  the  electrolytic  solution  on  the  current  density  v/ixh  fixed 
applied  overpotentiai.  From  the  three  figures,  we  can  see  thax 
the  current  density  is  approximately  proportional  to  xhe  bulk 
concentration  of  potassium  nitraxe,  except  for  the  case  of  -O’. 

Vfnen  electrolyxic  solution  ccnxains  cadmium  ions,  a 
differenx  sex  of  resulxs  was  obtained.  The  experimental 
currenx-density-vs . -time  curve  obtained  in  a  0,00251-1 
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Effect  of  tl\e  Bulk  Concentration  of  the  Electrolytic  Solution 
on  the  Current  Density  with  Fixed  Applied  Overpotential . 
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Cd(N0^)2  solution,  and  at  -G.40V  applied  overpotential  is  sho'/m. 
in  Figure  6.  In  addition  to  the  electrochemical  reduction  of 
the  nitrate  ions  on  the  electrode  surface,  a  precipitation 
reaction  also  takes  place  in  the  solution  because  the  cadmium 
ions  in  the  electrolytic  solution  can  coprecipitate  'vith  the 
hydroxy  ions  produced  by  the  electrochemical  reaction  on  the 
electrode  surface.  The  current  density  is  expected  to  be 
greater  than  that  when  the  cadmium  ion  is  absent  from  the 
electrolytic  solution.  Figure  10  shows  the  current-density- 
vs.-time  curves  obtained  in  a  O.OO5M  KT-iO^  solution  and  in  a 
O.OO25M  Cd(N0^)2  solution.  The  applied  overpotential  is  -O.i-OV 
for  both  cases.  From  the  curves,  one  can  see  that  the  current 
density  obtained  in  the  0.0025M  Cd(N0^)2  solution  is  larger 
than  that  obtained  in  a  O.OO5M  KNC^  solution.  The  reason  for 
this  is  that  the  consumption  of  the  hydroxy  ions  due  to  the 
precipitation  reaction  establishes  a  driving  force  to  produce 
more  hydroxy  ions  and  thus  promote  the  electrochemical  reaction. 
Some  interesting  facts  may  be  seen  from  the  figure.  Firstly, 
the  difference  between  two  current  densities  increases  as  time 
increases,  finally  reaching  a  constant  difference:  secondly, 
the  rate  of  decrease  of  the  current  density  with  time  is  slov/er 
for  the  case  of  using  Cd(N0^)2  as  the  electrolytic  solution. 

One  may  thus  obtain  a  larger  steady  state  current  density  using 
Cd(N0^)p  as  the  electrolytic  solution. 
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CHAPTER  5 

NUMERICAL  SOLUTIONS 
OF  THE  THEORETICAL  MODEL 

5 . 1  Orthogonal  Collocation  ?-Te~chod 

There  are  many  numerical  techniques  which  can  be  used 
to  solve  the  differential  equations  of  the  theoretical  model 
presented  earlier  in  Chapter  3* 

Orthogonal  collocation  method  was  selected  as  the  method 
of  solution.  The  orthogonal  collocation  method  requires  very 
little  implicit  mathematics  and  results  in  a  large  savings  in 
computer  time  over  the  common  finite  difference  scheme  such  as 
the  Crank-Nicholson  method.  A  brief  description  of  the  ortho¬ 
gonal  collocation  method  is  presented  in  Appendix  C. 

5*2  '/Vorking  Equations 


5*2.1  Dimensionless  Form  of  the  Differential  Equations 


The  orthogonal  collocation  method  requires  the  region  of 
solution  to  be  restricted  to  the  interval  C C , O  .  This  is 
accomplished  by  normalizing  the  spatial  variable  by  the 
parameter,^,  the  diffusion  thickness,  ^=x//. 


diniensionless  formas 

*  ■» 

P  —  - 

^  c"  ^  c° 

1  1 

the  dimensionless  time  is  defined  as: 


Substituting  the  above  dimensionless  variacles  into  equations 
(1?)  and  (19)  aJtd  rearranging  leads  to  the  following  dimension¬ 
less  eauations: 


^  * 

'"Z  ~^2  ^  ^  3?lst) 

or  '2'^"  ■  r  ‘ 


»hsre  K=  -i-£—  and  2Ksp=  , 

are  the  dimensionless  homogeneous  reaction  rate  constant  and 
solubility  product,  respectively.  Similarly,  equations  (21) 
and  (22)  become: 


^  ^2.  ♦ 


♦  ?  ♦ 


V>J 
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The  initial  conditions  are: 


1.  atr=0,  C.  =1, 


for  all  0 


if  i  1, 


2.  atr: 


=0,  €2"=  C2/CJ  ,  for  all  0  <  y  < 


1. 


The  first  two  boundary  conditions  come  from  the  dimensionless 
form  of  equations  (32)  and  (33) s 


1.  (-Vi-)  '  ^  ^ 


*  *, 


S  L‘2'>  J^=o 


K3  (c’“)3  (C°')3  -  (X.CJ)  cj’ 


=  0  , 


1 

2.  3(— ) 

C2  iTi)  J^=0 


r 


=  0 


1=0 


The  other  two  boundary  conditions  are: 


3 .  at  =1 ,  =1  , 


f  or  >  0  , 

4.  aty=l,  c^*=c^/cr.  for't>0. 


5.2.2  Discretized  Equations 


The  solution  is  approximated  by  a  (?.'+2 )  th-order  Legendre 
polynomial  which  satisfies  the  differential  equations  and 
boundary  conditions  exactly  at  N+2  points.  Those  coin's  are 
chosen  to  be  zeros  of  che  Legendre  polynomial  of  degree  T'  over 
the  interval  (0,1)  and  with  the  boundary  points  0  and  1. 
Discretizing  the  spatial  derivacives  ac  those  points  leads  to 
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the  following  set  of  firsx-order  ordinary  differential  equations 


i=  2, . N-rl,  (39) 


( iffQ ) 

where  ^  i  i=2,.,.  N  +  1  are  the  N  interior  collocation  points 

while  ,  and  exterior  collocaxion  points  which 

correspond  to  the  boundary  points  ax  =0  and  =1, 

respectively,  and  0=  '|  ^  <  i=2,  .  . .  N  +  1.  3^  . 

is  a  (N+2)  by  {N+2)  coefficient  matrix  which  comes  from  the 

discretization  of  the  second  derivative  {  or 

^  ^2  '  ^ '5  ^  each  collocation  poinx.  A  detailed  explanation 

of  matrix  3  is  sho'/m  in  Appendix  C.  G,  •I'f )  and  G,  ( %  , , 't ) 

+  ^  0  2  5  j 

are  xhe  concentrations  of  species  1  and  2,  respectively,  ax 
time  Z"  and  position 
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where  ^  is  aii  element  of  another  (N-^2)  by  (^n  +  Z)  coefficien' 
matrix  which  comes  from  the  discretization  of  the  first-order 
derivative  (  2)C^  or  )  at  each  collocation  points . 

Since  only  boundary  conditions  1  and  2  (equations  (43)  and  (44)) 
require  the  evaluation  of  the  first-order  derivative  at  ^ =0 , 
only  the  first  row  of  the  matrix  A,  i.e.,  A,  is  used  hers. 

^  f  o 

The  initial  conditions  for  both  cases  are  the  same: 

I.C.  's 

1.  at  -^=0,  C^*(^^,'^)=  1,  i  =  l,,,.N+2, 

2.  at  -^=0,  ^2/0^  ,  i=l,..,M+2, 


5 '2. 3  Calculation  of  Current  Densit\ 


The  current  density  is  related  to  the  flux  of  the  nitrate 
ion  at  the  electrode  surface  by  equation  (34).  The  discretized 
and  dimensionless  form  of  equation  (34)  is: 
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5 . 3  The  Introduction  of  a  Solir.e  Point:  ~o  tr.q  3i.;creti::ed 
Equations 


The  diffusion  thickness  introduced  in  the  derivation 
of  the  dimensionless  form  in  the  previous  section  can  he 
justified  hy  selecting  the  total  electrolysis  time  such 
that  diffusion  effects  are  negligible  at  x=i  during  that 
period.  At  the  beginning  of  the  electrolysis,  the  diffusion 
thickness  selected  according  to  the  total  electrolysis  time  is 
too  large.  This  leads  to  a  concentration  profile  which  is 
flat  in  most  of  the  interval  with  a  sharp  drop  in  a  very  small 
region  near  =0.  This  would  require  a  large  number  of 
collocation  points  leading  to  a  very  stiff  set  of  ordinary 
differential  equations. 

This  difficulty  can  be  overcome  by  the  use  of  a  method 

( 2.S) 

called  spline  collocation'  '  based  on  the  diffusion  boundary 
concept.  This  method  maintains  a  fixed  low  number  of  collo¬ 
cation  points  in  the  interval  v/hich  is  very  small  initially 
and  will  be  increased  as  time  increases  to  account  for  the 
thickening  of  the  diffusion  layer.  In  the  regions  outside 
this  interval,  the  concentration  is  flat.  As  a  matter  of 
fact,  the  spline  point  can  be  located  as  close  to  the  inter- 
face  as  to  achieve  any  desired  accuracy ' . 

t  the  concentration  gradient 


In  order  to  be  sure  tha 


at  the  spline  point  is  zero,  the  concentration  at  the  Nth 
collocation  point,  i.e.,  the  last  one  cefore  the  spline 
point,  is  tested  to  assure  that  it  is  within  a  very  small 
range  of  the  bulk  condition.  At  the  time  when  such  a  test 
fails,  the  spline  point  is  moved  further  into  the  solution. 

Introducing  a  spline  point,  0<  g  2.  new  variable, 

z='|’/ is  introduced.  The  discretized  equations  we  octained 
for  the  case  of  solution  containing  cadmium  ions  are: 


dC 

1 

df  12; 


h  lU  -  J 


i“2>  «  B 


±. 

=  (  >2  ) 

y  3;  -.C  (2;,T) 

-2K  [ 

2 .  T  s 

^  *  m  J.  f  J  C  J  J 

j  =  l 

1 

o  O 

2 


i=2 , . . .N*l . 


ror  the  case  of  solution  containing  no  cadmium,  ion, 


^  I 


Equations  (^4-3)  through  (46),  which  come  from  the  four  bcundar; 
conditions,  become: 

2  ^b  :i+2 

1,  (-|-  )  (  ^  izy^  ^=1'T  ) 

-  =0,  for  allt>0:  (i2) 


2.  3(' 


.  r-  ^ 


2  4  J-l  >  ,  C  J 


for  all  X>^* 


3.  ‘^i*'-:ii-2’t )'  -•  ail  X>o, 

•"•  h*<'|s'T>=  -2*<a:i,.z.t)=c|/3‘b  ior  m  t2C 


The  initial  conditions  become: 
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A 


T  r  f  . 

X  •  •  o  • 

1.  at  '^'=0,  iz^,^)=  1,  i  =  l,  .  .  .r;+2, 

2.  at  r=0,  C^*iz^,'^)=  C°/C°,  i  =  l,..,N+2. 


The  current  density  should  he  changed  to : 


5 Solution  Procedures 


For  the  case  of  a  solution  containing  cadmium  ions, 
equations  (48),  (49)  ar.d  (52)  through  (55)  i  with  the  initial 
conditions,  are  the  working  equations  used  to  solve  for  the 
concentration  profiles  of  NO^”  and  0K~  ions.  The  Cd"^^  ion 
concentration  is  calculated  by  the  electroneutrality  conditic 
in  the  solution: 
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For  the  case  of  a  solution  containing  no  cadmiun  ion, 
equations  (50)i  (51)  and  (52)  through  (55)  with  the  initial 
conditions,  are  the  working  equations. 

Equation  (56)  is  used  to  calculate  xhe  current  density 
at  each  time  both  cases. 

5 . 5  Computer  Program  Structure 


The  computer  program  for  solving  the  discretized 
equations  consists  of  a  main  program  and  5  rnain  subroutines. 
The  first  main  subroutine  evaluates  the  collocation  points 
(roots)  of  the  Jacobi  polynomial  of  order  N,  as  well  as  the 
firs  and  second  deri"atives  ^f  tne  polynomial  at  the  roots. 
These  derivatives  are  then  used  in  a  cecond  main  subroutine 
to  calculate  the  discretization  coefficient  matrices  A  and  3. 
The  third  main  subroutine  is  used  to  perform  a  semi-implicit 


integration  using  Gear's  routine 


(29) 


Tnis  suoroutine  caxis 


4  external  subroutines.  The  first  of  these  contains  the 

explicit  expressions  (48)  and  (49)  (or  (50)  and  (51)  in 

another  case)  for  the  discretized  coupled  first  order  differentia] 

equations  containing  the  3.  .  terms:  the  second  contains  the 

~  *  j 

Jacobian  matrix  for  the  same  equation:  the  third  performs  the 
first  stage  of  Gaussian  eliniinatior.  '  of  Jacobian  matrix; 
tne  _ast  performs  the  back  substitut.cn  cf  laussiar.  eiimi..ation 


of 
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of  the  Jacobian  natrix.  The  first  external  subroutine  called 
by  the  integration  subroutine  also  calls  another  subroutine 
coming  from  the  IMSL  package.  This  subroutine  solves  equations 

(52)  and  (53)  simultaneously  to  give  the  surface  concentrations 

*  * 

of  C.  and  C_  .  The  IMSL  subroutine  calls  another  function 

1  c 

which  feeds  the  equations  to  be  solved.  Before  writing  the 

* 

function,  (surface  concentration  of  )  is  solved 

in  terms  of  i  =  l ,  .  .  ,  and  (t-i'lT)*  i  =  l ,  .  .  , 

using  equation  (53)»  and  then  substitute  into  equation  (52) 

♦ 

to  eliminate  C2  Thus,  only  one  equation  which 

■» 

contains  only  one  unknown,  is  needed*  Cnee  v/e 

obtain  (z^,'^),  can  be  calculated  easily. 

The  Jacobian  matrices  for  both  cases  can  also  be  evalu- 

■» 

ated  after  substituting  the  expression  of  into 

equations  (43)  and  (49)  (or  (50)  and  (51)  in  another  case). 

For  the  case  of  solution  containing  no  cadmium  ion,  the  Jacobian 
matrix  is: 
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which  is  a  2M  by  2N  .~atrix.  For  the  case  of  solution 

cadmium  ions,  the  Jacobian  matrix  is  the  same  except 

the  terms  in  the  diagonal  line  of  the  third  quadrant 

subtracted  by  K  and  the  terms  in  the  diagonal  line  of 

* , 

fourth  quadrant  are  all  subracted  by  ZKif^iZDKs'p/C^ 


containing 

that 


are  all 
the 


r) 


The  values  obtained  from  the  integration  subroutine  are 
then  used  in  a  fourth  main  subroutine  to  calculate  the  current 
density  by  equation  (56).  The  last  main  subroutine  is  used  to 
evaluate  the  concentration  at  any  desired  point  betv/een  0  and  1 
by  interpolation  using  the  values  at  the  collocation  points. 

The  structure  of  the  whole  program  is  shown  in  Figure  11. 
Programs  from  both  cases  are  listed  in  Appendix  D,  Programs 
herein  were  executed  on  an  AIOAHL  470/V6  computer. 
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fitti?:g  of  the  theoretical 

MODEL  WITH  EXPERIMENTAL  DATA 

6 . 1  The  PeterTiiination  of  the  Electrochemical  .-.eaction 
Kinetics  ?ar3-’ne~ers 

The  unkno’-vn  heterogeneous  reaction  rate  cons¬ 
tants  which  appear  in  equation  (52)  were  determined 
from  the  current-density-vs . -time  data  obtained  in  a 
solution  which  did  not  contain  cadmium  ion.  The  di¬ 
fferential  equations  (50)  and  (51)  which  describe  the 
transport  process  were  solved  by  computer  to  determine 
the  rate  constants.  Comparing  equation  (3i)  with  equa 
tion  (31a),  one  can  see  that  both  K^  and  K,  depend  on 
the  applied  potential  with  exponential  dependence. 
Thus,  X.  and  K2  should  be  constant  when  the  applied 
potential  is  maintained  at  a  constant  value  during  the 
electrolysis.  K,  and  X^  were  obtained  by  comparing  the 
experimental  data  in  the  absence  of  cadmium  ion  in  the 
electrolyte  with  the  computer-calculated  current  densi 

The  experimental  current-density-vs , -time  data 


data 


with  the  electrolyte  containing  0.005'^  nitrate  ions 
and  the  applied  overpotential  of  -0.40V  was  first 
fitted  to  the  computer-siaulated  data  to  obtain  the 
and  values.  Due  to  the  lack  of  information 
about  the  values  of  and  K21  rough  estimates  of 
and  K,  were  evaluated  by  substituting  the  applied 
potential  V,  which  is  equal  to  the  value  of  the  applied 
overpotential  plus  the  rest  potential,  to  equations 
(31)  (jla).  The  estimated  values  of  X-  and  K,  were 

increased  or  decreased  until  the  computer-calculated 
current-density-vs . -time  curve  had  the  same  shape  as 
the  same  curve  obtained  experimentally,  (hote,  the 
background  water  decomposition  current  has  been  sub¬ 
tracted.  ) 

After  we  obtained  the  right  shape,  the  anodic 
reaction  rate  constant,  K^,  was  then  fixed.  5y  changing 
K, ,  one  can  make  the  entire  curve  shift  upward.  Con¬ 
versely,  decreasing  the  X^  value  will  decrease  the 
current  density  at  each  time  point  and  thus  shift  the 
entire  curve  downward.  A  set  of  current-density-vs . - 
time  curves  with  fixed  X^  values  and  different  values 
of  X.  are  plotted  in  Figure  12.  The  values  of  X,  and 
X,  when  the  concentration  of  nitrate  ions  is  O.CCfVi 
and  the  applied  overpotential  is  -G,-»GV  v/as  obtained 


Eleutrolyiiis  Time,  mt; 
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by  interpolating  the  experi.-en'cal  curve  cecv/esn  chose 
curves  in  Figure  12.  The  values  of  and  '-vhich 
achieved  the  best  approximation  of  the  experimental 
data  aire  K^=  0.12  ^10~^  and  K2=  0.6  5<10~-^,  respectively 

These  values  were  chen  used  co  predicc  the  cu- 
rrent-density-vs . -time  curves  at  che  ocher  ion  con¬ 
centrations.  Figure  13  shows  the  computer-simulated 
current  density  curves  and  experimental  curves  ac  va¬ 
rious  nitrate  ion  concentracions . 

The  values  of  and  at  the  applied  overpo¬ 
tential  of  -0.60V  were  different  from  those  ac  -0.4CV. 
The  sane  procedure  was  used  in  evaluating  the  values  of 
and  at  the  condition  of  -O.oOV.  Figure  Ilf-  shows  a 
set  of  calculated  current-density-vs . -time  curves  with 
fixed  valued  and  different  values  of  .  The  value  o 
at  the  condition  of  -0,60V  was  then  obtained  in  the 

i. 

same  way  as  before.  The  values  of  and  K2  were  deter- 
mined  to  be  0.10  x  10  and  C.50X.IC  ,  respectively. 
These  values  were  then  used  to  predict  the  current-den¬ 
sity-vs  .  -tine  curves  at  the  other  ion  concentrations . 
Fig’ire  15  shows  the  comparison  of  the  predicted  ar.d  ex¬ 
perimental  current-density-vs . -time  curves. 

The  reaction  rate  constants  ac  the  applied  overpo 
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Fi^Mue  I'Jt  Exi)ei-imeiital  and  Computei'-Simula ted  Curren t-Densi ty-vs , -Time 

Curves  with  Varioua  NiLi-ato  Ion  Bulk  Coticentrationa  in  the  Caae 
of  r:i  fi;  I  I'ol  yte  CUjntainiri'f  No  Cadiiiium  Ion. 
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tential  of  -O.SOV  were  not  pursued.  The  reason, which 
was  discussed  previously,  is  that  az  this  high  applied 
overpotential  iriore  complicate  sequences  of  reactions 
take  place.  This  simple  reaction  expression  can  not 
describe  the  phenomena. 

The  heterogeneous  rate  expression  so  determined 
is  believed  to  be  a  correct  one.  This  claim  is  supported 
by  the  fact  that  the  rate  constants,  and  K2,  are 
independent  of  the  bulk  concentration  at  a  given  poten¬ 
tial.  This  is  indeed  the  case  as  v/as  shewn  in  Figures 
13  and  15. 

6.2  The  Determination  of  the  Homogeneous  Frecititation 
Reaction  Rate  Constant 

Higher  current  density  was  observed  when  cadmium 
nitrate  was  used  instead  of  potassium  nitrate  as  tne 
electrolyte  in  the  electrochemical  cell.  The  hydroxy 
ion  produced  by  the  reduction  of  nitrate  ict.  copreci¬ 
pitated  with  the  cadmium  ion  in  the  solution.  This 
increased  the  reduction  reaction  rate  and  led  to  higher 
current.  This  homogeneous  reaction  was  assumed  to  be 
linearly  dependent  on  the  degree  of  supersaturation  of 
cadmium  hydroxide  (see  Squation  (1^)}.  Since  the  heza- 


roger.eous  reactior.  rate  constants  K.  and  K-,  v;ere  deter- 
~ined,  the  homogeneous  rate  constant  d  could  be  evalu¬ 
ated  'ey  fitting  the  second  set  of  experixtental  data 
with  the  theoretical  model  to  determine  tne  value  of  1. 

Fig’ure  16  shows  the  current-density-vs , -time 
curves  with  various  jc  values  while  the  other  parameters 
were  held  at  a  constant  value.  The  current-density-vs . - 
time  data  o'etained  in  a  0.002j7Ji  cadmium  nitrate  solution 
at  the  applied  overpotential  of  -O.^OV  is  also  plotted  in 
rig’ure  16.  The  reaction  rate  constant,  >,  was  then  eva¬ 
luated  approximately  from  Figure  16,  The  most  pro'oa'cle 
homogeneous  reaction  rate  constant  was  determined  to  'oe 
50  sec“^. 


k=  100 


ity-vti . -Time  Curveii  with  Fixed 
ut;  Reaetioii  Rate  ConatairLs. 
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CHAPTER  7 

DISCUSSIONS  AT^D  CONCLUSIONS 

7 •  1  Surface  Concentrations  of  Various  Ions  5.s  Runctions 
of  Time 


The  kno'.vledge  of  the  surface  ion  concentrations  can  pro 
vide  information  concerning  the  electrochemical  deposition 
process.  Although  the  surface  ion  concentrations  -.vere  not 
available  from  the  experiment,  that  information  could  be 
generated  by  computer  simulaxions. 

Figure  17  shows  the  surface  concentrations  of  NO^"  and 
0H~  ions  as  functions  of  time  in  the  absence  of  cadmium  ion 
in  the  solution  at  t'wo  potentials,  -O.LO’/  and  -0.60V  from 
equilibrium  potential  and  at  a  N0^~  bulk  concentration  of 
0.CO5M.  In  general, the  surface  concentration  of  N0^~  ions 
decreases  as  time  increases,  while  the  surface  concentration 
of  0H~  ions  increases  as  time  increases.  This  is  due  to  the 
electrochemical  reaction  (23),  in  which  one  nitrate  ion 
reacts  with  two  electrons  and  which  produces  three  hydroxy 
ions . 

For  tne  case  of  -O.m-CV  overpotential,  the  CH”  concen¬ 
tration  increases  sharply  from  the  bulk  value  of  nearly  cere 


KEY 


74 


75 


to  abouo  0.5  dimensionless  concenxraticn  unins  in  less  nhan  1 
msec  before  the  hydroxy  ions  generated  are  removed  effecniveiy 
by  diffusion.  Thereafter,  the  0H~  concentration  increases  at 
a  slower  rate.  For  ions,  the  surface  concentration  drops 

to  about  0.72  dimensionless  concentration  units  from  the  bulk; 
value  of  1.0  in  less  than  1  msec  and  then  decreases  at  a 
slower  rate  when  the  diffusion  can  effectively  supply  the 
reactant  for  the  electrode  reaction.  The  sharp  char.ges  of  zhe 
surface  concentrations  of  NO^  and  CH  icns  cause  the  snarp 
decrease  of  the  current  density  in  a  short  time  period  in 
the  beginning  of  mhe  electrolysis. 


At  a  more  negative  overpotential  of  ~C.60V,  the  increase 
in  surface  CH”  concentration  and  the  decrease  of  NO^  concen¬ 
tration  are  even  more  significant  and  last  over  a  longer  time 
period.  The  surface  concentration  of  drops  to  about 

0.5  dimensionless  concentration  units  and  CH  concenxration 
increases  to  about  0.S7  dimensionless  concentration  units  in 
2.0  msecs.  The  behavior  is  xhe  consequence  of  the  higher 
overall  reaction  rate  at  -O.6OV.  The  nitrate  ions  were  con¬ 
sumed  at  a  higher  rate  which  produced  more  hydroxy  ions, 
longer  induction  time  is  needed  for  diffusion  zo  effecxively 
supply  the  reactant  from  the  bulk  soluxion  and  remove  xne 
oroducx  from  the  elacxrode  surface. 


Similar  profiles  of 


suriact 


cor 


and  CH 
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at  lower  bulk  N0^~  concentra~ions  were  also  obtained.  These 
are  shown  in  Figures  18  and  19* 

■;ifhen  the  electrolyte  contains  cadniun  ions,  the  homo¬ 
geneous  precipitation  reaction  of  cadmium  ions  and  hydroxy 
ions  near  the  electrode  surface  promotes  the  electrochemical 
reaction  on  the  electrode  surface  in  the  cathodic  direction. 
This  causes  the  nitrate  ions  to  be  consumed  at  a  higher  rate 
than  the  case  when  no  cadmium  is  present.  This  is  shewn  in 
Figure  20  where  the  homogeneous  reaction  rate  constant  is  50 

4 

sec"'^.  It  is  interesting  to  note  that  the  cadmium  ion  surfac 
concentration  increases  initially  and  then  starts  to  decrease 
This  phenomenon  is  explained  as  follows.  Immediately  after 
the  electrolysis  began,  ions  were  consumed  by  the  elec¬ 

trochemical  reaction  which  produced  an  NO^"  concentration 

gradient  made  the  NO^'  ions  in  the  bulk  solution  diffused 

+2 

toward  the  electrode  surface.  The  Gd  ions,  on  the  other 

hand,  moved  with  the  NO^'  ions  in  the  same  direction  in  order 

to  maintain  the  elecroneutrality  condition.  In  the  meantime, 

the  OH  concentration  was,  however,  not  so  high  as  to  consume 
+2 

all  the  Cd  ions  brought  in  by  the  transport  process.  The 
cadmium  ions  was  thus  accumulaued  arid  resulted  in  a  higher 
value  than  its  bulk  condition.  As  scon  as  the  surface  concen 
tration  of  OH  ion  raised  to  some  value,  tne  surface  concen¬ 
tration  of  the  Cd  ^  iens  starts  to  decrease  because  of  the 
abundant  supply  of  OH"  ions  consuming  the  cadmium  ions  at 
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this  time. 


For  the  nitrate  ions,  the  surface  concentration  also 
decreases  as  time  increases  at  a  rate  somewhat  larger  than 
that  for  the  case  of  solution  containing  no  cadmium  ion.  This 
is  shown  in  Figure  21.  Figure  21  also  shows  the  difference  in 
the  surface  concentration  of  the  hydroxy  ion  between  those  two 
cases,  i.e.,  in  the  presence  and  in  the  absence  of  cadmium 
ions.  The  initial  behaviors  of  two  cases  are  very  similar. 
After  this  induction  period,  the  ion  concentration  for  the 
case  of  a  solution  which  contains  cadmium  ions  decreases  at 
a  more  gradual  rate  as  a  result  of  the  homogeneous  precipi¬ 
tation  reaction  between  cadmium  and  hydroxy  ions. 

7.2  Concentration  Profiles  of  Various  Ions  at  Various  Time 


For  the  case  of  solution  containing  no  cadmium  ion,  the 
nitrate  ions  were  consumed  and  the  hydroxy  ions  were  produced 
on  the  electrode,  the  concentration  gradients  were  established 
which  made  the  nitrate  ions  diffused  from  the  built  solution 
toward  the  electrode  surface  and  the  hydroxy  ions  diffused 
away  from  the  electrode  surface  toward  the  built  solution. 
Figures  22  through  27  shov/  the  concentration  profiles  of  HC:," 
and  OH  ions  with  different  applied  overpotentials  at  some 
selected  times.  From  these  figures,  one  cari  see  that  tne 
diffusion  thicicness,  ^  ,  increases  v/imh  time.  The  diffusion 
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I'i^rure  22  Goacenti’ation  Profiles  of  Nitiate  and  Hydroxy  Ions  at 
Various  Time  for  the  Case  of  Solution  Containing  No 
Cadmium  Ion. 
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Cadiiiiurii  Ion. 
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Fi  trure  26 1  Concen  tration  Profiles  of  Nitrate  and  Hydroxy  Ions  at 
Various  Times  for  the  Case  of  Solution  Containing  No 
Cadmium  Ion. 
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centration  Pfofiles  of  Nitrate  and  Hydroxy  Ions  at 
ions  Times  foi'  the  Case  of  Solution  Containing 


thicknesses  at  various  tine  are  listed  in  Tables  2  and  3> 

The  concentration  profiles  of  N0^~,  OH”  and  ions 

at  various  time  for  the  case  of  solution  containing  cadmium 
ions  are  shown  in  Figures  28  and  29.  The  shape  of 
concentration  profiles  is  the  same  as  that  for  the  case  of 
containing  no  cadmium  ion  but  the  ?10^~  concentration  on  "che 
surface  is  somewhat  lower  than  that  for  the  case  of  contain¬ 
ing  no  cadmium  ion.  The  diffusion  thickness  increases  as  time 
increases . 

In  the  previous  discussion  of  the  surface  concentration, 
the  surface  concentration  of  OH”  ions  increases  in  the  beginn¬ 
ing.  During  that  time,  OK”  ions  diffuse  from  the  electrode 
surface  to  the  bulk  solution.  Figure  28  shov/s  that  the 
surface  concentration  increases  and  the  diffusion  layer  thick¬ 
ness  also  increases  as  time  increases.  Then  the  surface 
concentration  starts  to  drop  and  the  diffusion  layer  thickness 
stays  at  about  the  same  value.  This  is  sho'wn  in  Figure  30. 
Table  4  shows  the  diffusion  layer  thickness  for  the  diffusion 
of  OH”  ions  as  a  function  of  time.  In  Table  4,  one  can 
see  that  the  diffusion  layer  thickr.ess  stops  increasing  as  a 
result  of  the  decreasing  surface  concentration  of  OK”  ion. 

+2 

The  Cd  concentration  profiles  are  snov/r.  in  Figure  29  > 
As  the  profiles  indicate,  for  each  specific  time,  when  the 
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Table  2:  Diffusion  Thickness  as  a  Function  of  Time 
for  the  Case  of  Solution  Containing  'Jo 
Cadmium  Ion.  The  Applied  Cverpoten'ial  is 
-0.40V.  The  Nitrate-Ion  Bulk  Concentration 
is  0.005M. 


Time  (sec) 


0.1235  Xio 


-5 


0.6121  X 10 


-5 


0.1758  XIO 


-4 


0.3572  XIO 


-4 


0.1033  XlO"^ 


0.4737  SIO"^ 


0.1258  XlO 


-2 


0.7417  xio' 


0.3946  XIO 


-1 


0.7032  XIO 


-1 


Diffusion  Thickness( cm) 


0.25  XIO 


-4 


0.45 


0.80  XIO 


—  •4' 


0.15  X 10 


0.25  X 10 


-3 


0.50  X 10"^ 


0.80  XIO' 


0,20  XIO 


1 


0.45  X 10 


-2 


0,60  X 10 


-2 
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Table  3:  Diffusion  Thickness  as  a  Function  of  Time 
for  the  Case  of  Solution  Containing  I'fo 
Cadmium  Ion.  The  Applied  Overpoten'cial  Is 
-0.60V.  The  Nitrate-Ion  Bulk  Concentration 
Is  0.005M. 


Time  (sec) 


0.1313  xio"-5 


0.4982  XlO 


-b 


0.1642  XIO 


0.3318  X10‘ 


0.1298  xio"^ 


0.4407  XIO”^ 


0.1258  XIO 


-2 


0.8718  XIO 


-2 


0.3229  XIO 


-1 


0.7004  XIO' 


Diffusion  Thickness  (cm) 


0.25X  10 


-4 


0.40  X 10 


0.70  X  lo""*’ 


0.10  X 10' 


0.25  XIO 


0.45  xio"^ 


0.75  xio"^ 


0.25  XlO 


-2 


0.45  XIO' 


0.65  XIO 


-2 
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30 «  Concentration  Prol  ilea  ol’  Hydroxy  Ion  When  the  Surface 
Concentration  of  Hydroxy  Ion  Starts  to  Dxop. 
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Table  4:  Diffusion  Thickness  of  OK”  Ion  as  a  Function 
of  Time  for  the  Case  of  Solution  Containing 
Cadmium  Ions.  The  Applied  Overpotential  Is 
-0.40V.  The  Nitrate-Ion  Bulk  Concentration  Is 
0.005M.  The  Cadmium-Ion  Bulk  Concentration  Is 
0.0025M. 


Time  (sec) 

Diffusion  Thickness ( cm) 

0,1236  xio”-^ 

0.25  XIO”*^ 

0.6128  XIO"^ 

0.45  xio”^ 

0.1798  xio"^ 

0.80  XIO”^ 

0.3513  xio’^ 

0.15  XIO”^ 

0.1017  X 10"^ 

0.25  XlO"^ 

0.4856  X 10"^ 

0,50  X 10”^ 

0.1148  XIO”^ 

0.70  xio"^ 

^  0.8811  XIO"*^ 

0.85  xio”^ 

0.4050  X 10”^ 

0.85  xio”^ 

0.7031  Xio"^ 

0.35  X 10”^ 
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concentration  near  the  electrode  is  greater  than  the  bulk 

■*•2 

concentration  of  Gd  ,  its  value  versus  distance  from  the 
electrode  decreases  slowly  in  the  beginning,  then  faster,  and 
more  slowly  again  before  finally  reaching  its  bulk  value.  At 
the  time  when  the  concentration  values  near  the  electrode 
surface  are  stili  lar-ger  than  the  bulk  concentration,  the 
concentration  decreases  to  some  minimum  value  at  some  distance 
from  the  electrode  surface  and  then  increases  until  it  reaches 
the  bulk  value.  The  location  of  the  minimum  concentration 
point  moves  toward  the  electrode  surface  as  time  increases. 

At  sufficiently  long  time,  t.he  minimum  point  finally  reaches 
the  electrode  surface. 

7 . 3  Conclusions 


The  transport  model  for  the  electrode  kinetics  and  the 
homogeneous  reaction  presented  previously  has  already  explained 
very  successfully  the  characterestics  of  the  current  density  as 
a  function  of  time  for  both  the  case  of  solution  containing 
cadmium  ion  and  that  without  cadmium  ion.  This  model  was  then 
used  to  simulate  the  concentration  profiles  for  various  ions 
involved  in  the  reactions.  The  behavior  of  the  various  ions  on 
the  electrode  surface  and  in  the  solution  was  obtained.  The  di 
fusion  processes  v/ere  also  knovm  after  examining  these  concen¬ 
tration  profiles  carefully.  These  can  be  summarised  as  follov/s 


1»  The  current:  density  decays  with  time  aiter  a  double-layer 
charging#  The  time  for  this  charging  is  so  short  that  v/e 
can  not  observe  it  by  the  experiment.  The  decay  for  the 
current  density  is  a  result  of  the  depletion  of  the  nitrate 
ions  on  the  electrode  surface. 

2.  The  current  density  for  the  case  of  a  solution  containing 
cadmium  ions  is  greater  than  that  for  the  case  of  solution 
containing  no  cadmium  ion.  This  is  because  the  precipitaticr 
of  Cd(CH)2  promotes  the  electrochemical  reaction  cn  the  elec" 
rode  surface  by  consuming  0K~  generated  in  this  reaction. 

3.  The  current  density  depends  on  the  bulk  concentration  of 


nitrate  ions  and  the  applied  overpotential.  Increasing 
the  NO^"  bulk  concentration  or  the  applied  overpotential 
will  increase  the  current  density. 

The  applied  overpotential  is  one  of  the  factors  which 
changes  the  reaction  rate  of  the  electrochemical  reaction 
on  the  electrode  surface.  Increasing  the  applied  over¬ 
potential  will  increase  the  overall  reaction  rate. 

5.  For  the  case  of  solution  containing  no  cadmium  ion,  the 

surface  concentration  of  decreases  as  time  increases 

while  the  surface  concentration  of  OH""  increases  as  time 
increases.  This  is  because  ions  are  consumed  to  produce 

the  oh”  ions  during  the  electrochemical  reaction  on  the 
electrode  surface. 

6,  For  the  case  of  solution  containing  cadmium  ion,  the 

surface  concentration  of  icn  also  decreases.  The 

decreasing  rate  is,  however,  faster  than  that  for  the  c-'- 
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of  soluij.i..''.  containing  no  iad.T.iu.Ti  ion.  Tiie  surface 
concentration  of  OH  ion  increases  in  the  beginning  ar.d 
reaches  a  maximum.  When  the  production  rate  of  CH~  due  to 
the  electrochemical  reaction  is  less  than  the  consumption 
rate  of  0K~  due  to  the  precipitation  reaction,  the  surface 
concentration  of  CK~  starts  to  decrease.  To  maintain  elec¬ 
troneutrality,  the  surface  concentration  of  Gd"^^  ions  v/ill 
exceed  its  bulk  value  and  then  decrease  to  the  value  lower 
than  its  bulk  value. 

7.  The  diffusion  processes  for  the  case  of  solution  containing 
no  cadmium  ion  is  that  the  nitrate  ions  diffuse  from  tne 
bulk  solution  toward  the  electrode  surface  while  the 
hydroxy  ions  diffuse  in  the  opposite  direction. 

3.  For  the  case  of  solution  containing  cadmium  ions,  the 

diffusion  directions  for  both  the  nitrate  ions  and  hydroxy 
ions  are  the  same  as  that  for  the  case  of  solution  containing 
no  cadmium  ion.  However,  the  cadmium  ions  diffuse  toward 
the  bulk  solution  in  the  beginning  and  then  finally  char.ge 
direction  toward  the  electrode  surface. 

7  r.e  commend  at  ion  for  Future  .'/orks 

Although  the  model  has  successfully  explained  the  cases 
of  -0.GOV  and  -O.oOV,  is  is  not  likely  valid  for  the  case  cf 
-O.tCV  due  to  the  unexpected  characterestics  cf  tne  current- 
density-vs . -time  curves  and  tne  comparison  of  the  current 


density  '.vith  the  limiting  case.  V/hen  the  applied  cverpoter.tial 
is  as  high  as  -0.80V,  the  reaction  mechanisms  and  electrode 
behavior  should  be  identified  by  doing  some  further  experimencs 
Furthermore,  the  current-density-vs . -time  data  at  -C.cOV  applie 
overpotential  for  the  case  of  solution  containing  cadmium  ions 
is  not  available  yet.  It  can  be  obtained  by  repeating  the 
experimental  work  we  decribed  before.  Once  the  data  is  obtained 
the  model  can  be  used  for  comparison  with  the  experimental 
data. 


Finally,  the  configuration  of  the  electrode  in  the 
impregnation  process  is  not  as  simple  as  we  have  used  in  derivi 
the  model.  It  is,  in  fact,  a  porous  electrode  in  the  real 
case.  If  a  micro-porous  electrode  is  possible  to  obtain, 
the  experiment  can  be  repeated  to  obtain  the  current- density- 
vs.-time  data.  Cnee  the  data  are  available,  the  transport 
processes  can  be  examined  for  the  condition  of  porous  electrode 
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Contonts  lnst.rii(^t.ions  Cominerit ; 


oy'i)  C9  RCT  Heturti  to  main  progi'aiii 


APPENDIX  3 

LISTING  OF  ZXPERIMET.TAL  DATA 
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:a3i.a 


5.3  Currant-Der.si-y-vs  .-Tir.e  Za-a  for  the  Case  of 
Solution  Containing  '<'o  Cadniuti  Ion.  .'-'itrata- 
lon  3ulk  Concentration=  O.OCtM,  Initial  prl=  3*2, 
Aopliad  Overpotential=  -O.SOV,  Rest  Potential 
='  -0.35^V. 


01:3 
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Table  3. 


^  urr  ar.b  - 1  er.  s  i  ~y  -  vs 


SoluTiion  Gonbair.ir.g 


-ir.e  lor  "ne  ^a 

0  Cad.T.iuni  Ion.  .'iitr 
Ion  Bulk  Conc9ntra“ior.=  0.C025M,  Initial  pH=3.C, 
Applied  0v9rnotential=“0 . 6CV,  Rest  Potential 

=  -0.35^’^* 


(iJ 
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(\)  P  Pi  P 
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:atle  5.9  Curr9nt-5er.sity-ys,-ri.T.e  5a-a  for  -r.e  Tase  of 
Solution  Conlainir.^  r.'o  Cad.T.iu~  lor.,  iliorao 
Ion  Bulk  Concentra-:ion=  C,CC125M,  I.nioial 
Applied  0vert3o-e.n-;ial=  -O.cOV,  ?.eso  Pooenoia- 
=  -0.35^V. 


CUP F'Ef'^T  I'E’r ^5 1 T' 
'  nf'lF'  C  ‘’VZ 


.  000  1 
.  0002 
.  0003 
.  0004 
.  0005 
.  0006 


.  000  >■ 


.  0000 
.  000 '5 


1S4. 043 
141. 976 
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1 1 1 . 1 77 
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.012 
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.016 

13.4043 

.  0 1 3 

17.2776 

.  02 

16.5263 

•  ^  ^  7*  1  ” 

1 


••cJ  lli 
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Table  3.12  Background  Current  Density  of  Water.  Initial 
pK=  3.0,  Applied  Overpotential=  -0,307,  Rest 
Potential=  -0.35^’^* 


.  0007 

47. 7013 

.  0003 

43.9453 
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41.313 

.  00 1 

33.3333 
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37 . 1 344 
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35 .  ►^■32 

.  00 1 3 

34 . 1 793 

.  00 1 4 

33.4234 
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32.3013 

.  00 1 6 

31.5504 

.  00 1 7 

30 . 7992 

.  00 1 3 

30 . 043 

.  00 1 9 

29. 3724 

.  002 

23.9212 

.004 

22.533 

.  003 

19. 9033 

.  003 

13.0233 
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13.5233 

.012 
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.014 

14.3433 

.013 

13.3972 

.013 
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.  02 
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12. 3947 
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11.  ':'433 

.  023 

1 1 . 3433 

.  023 
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.  03 

1 0 . 3923 

.  0-^2 

10.5137 

.034 

10.5137 

.  033 

9. 73531 

.  033 

9. 73531 

.04 

9.  39 

cl  O 


123 


Tacle  3.16  Curren“-Der.sity-’/: 


;  o: 


Lse  01 


Solu-ion  Contair.ir.g  Cad.T.iun  Ions.  .'linraoe- 
lon  Buil-c  Cone  en  nr  a*;!  on  =  0 .0025Mi  _  Cadmi 
Sulk  Concennranicn=  0.00125M>  Initial 

Apiolied  Overpotential=  -O.^CV,  Has*  ro - 

=  -0.35^V. 


j  p.  p 


O 
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Tacle  5.19  Currer.t-Density-vs . -T 

Solutior.  Containing  Cadniun 


ata  for  “ne  Case  ci 
ms.  Mi^rate- 
lon  Bulk  Concen~nation=  0iCC125‘*ii  Cadniiun-Icr. 
Bulk  Concentration^  O.COO625M,  Ini~aal  pH=  3*0» 
Apnlied  Over‘ooten'Cial=  -Oi^OV,  Rest  r'O'jential 
=  -0.35^V. 


>--1 


Bulk  Cone sn'cra‘Sion=  0.00625M,  Initial  pi--- 
Auplied  Over'no'!:ential=  -0  ,Q0'f ,  Best  Potential 

=‘-Q.35^^* 


i-i 
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The  method  of  weighted  residuals,  ,  is  a  general 
method  for  obtaining  solutions  to  equations  of  change,  in 
our  case.  Pick's  second  law.  In  the  MWR ,  one  assumes  a  trial 
function,  usually  a  set  of  weighted  polynomials , substitutes 
this  trial  function  into  the  differential  equation  and  then 
selects  the  coefficients  of  the  polynomial  terms  by  specify¬ 
ing  that  the  residual  be  zero,  on  the  average,  az  certain 
points.  If  one  evaluates  the  differential  equation  at  the 
zeros  of  an  orthogonal  pol^'momial,  the  residual  will  of 
necessity  be  exactly  zero  at  these  collocation  points.  3y 
increasing  the  number  of  collocation  points,  the  trial 
function  would  satisfy  the  differential  equation  at  mors  and 
more  points^  ^  ^ , 


(1)  3. A.  Pinlayson  i  Richard  3ellman  (Ed.),  'Mathematics 

in  Science  and  Engineering',  Vol.  8?,  Academic  Press, 
New  York,  1972,  pp  9. 
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v/ith  the  bouTidary  conditions  of  -he  forr. 

y(0,t)=g(t)  (3) 

y(l,t)=h(t)  (4) 

and  initial  condition  such  as 

y(x,0)=p(x)  (5) 

For  a  mass  transfer  problem  D  would  be  a  dimension¬ 
less  diffusion  coefficient  and  f(y)  would  be  homogeneous 
reaction  term<  Geometric  considerations  indicate  that 

oC-^=0  and  mathematical  considerations  indicate  that  a 

( 2 ) 

suitable  but  not  unique  trial  function  is^ 

N 

y(x,t)  =  (l-x)y(0,t)+xy(l,t)+x(l-x)2Za  .(t)?,  .(x)  (6) 

■T_  I 

' 

The  ? j  ,  of  equation  (6)  refers  to  one  member  of  a 
complete  set  of  polynomials.  Table  C,1  lists  some  Legendr 
pol;ynomials ,  the  special  case  of  equation  (1)  v/here  ot  =  3 
and  their  roots  .  Note  that  a  Legendre  polynomial  of  iegr 
has  N  real,  distinct  roots  or.  the  interval  0<x<l, 


1.2)  E.k.  Finlayson  -  F.ichard  Fellm.a^-,  '  Zd 
in  Science  and  Zngineerin.g‘ ,  Vo*.  37. 
New  York,  1971,  pp  105- 


' Mathematics 
cdemiic  Press, 
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TA3L2  C.  1 


Legendre  Polynomial  and  Their 


:^.oots 


(3) 


n 

^N-l 

Roots  of  P., 

iN  ' 

1 

-l+2x 

0.500000000 

2 

l-6x+6x^ 

0.211324865 

0.788675135 

3 

-l+12x-30x^+20x^ 

0.112701665 

O.5OCOOOOOO 

0.887298335 

4 

l-20x+90x^-l40x^+70x^ 

0.069431844 

0.330009478 

0.669990522 

0.930568156 

(3)  vTohn  Villadsen,  'Selected  Approximation  Method  for 
Chemical  Engineering  Problems',  Heproset,  Copenhagen, 
Denmark,  1970,  pp  A3  and  A4. 
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I~  is  evident  that  regardless  of  the  value  of  M  cho¬ 
sen  one  can  always  rewrite  eqn.(6)  in  the  general  form 

M+2 

y(x,t)»y’  x^"^b..(t)  (7) 

<j  =  l 

One  could  at  this  point,  as  is  normally  the  case  in 
using  MWR,  obtain  the  coefficients  b.(t)  by  substitution  of 

j 

eqn.(7)  in  the  differential  equation  whose  solution  is  de¬ 
sired.  However,  the  computer  programming  is  greatly  simpli¬ 
fied  when  they  are  written  in  terms  of  the  solution  of  the 
differential  equation  at  the  collocation  points,  y(x4,t) 
where  i*l ,  2 , 3  •  •  tH+l ,  and  each  is  one  of  the  roots  of 
the  polynomial  or  one  of  the  two  boundary  values.  So  in 
terms  of  the  ?i+2  collocation  points  we  may  write.- 

y(x,  ,t)=  VxP'^b^vt)  (8) 

xi 

'.Ve  now  have  an  approximating  polynomial  function  which 
v/e  will  force  to  fit  our  partial  differential  eqn.(2)  at 
the  N+2  collocation  points  (i.e.  cnoose  the  b4(t)  so  that 
our  partial  differential  equation  is  satisfied).  In  order 
to  do  this  we  will  need  both  the  first  and  second  derivatives 
witn  respect  to  x.  These  derivatives  can  be  calculated  ex¬ 
plicitly  in  terms  of  x  since  we  have  already  defined  the 
dependence  of  y  upon  the  spatial  coordinate , eqn .( o ) and  (5). 

It  is  easily  seen  that  if  one  considers  all  the  co- 


13^ 


llocation  points,  equation  (5)  represen's  one  n'.enicer  o 

a  set  of  equations: 


N+2 

y(x^ , t)-  y~xj~^b j( t)  (9) 

0=' 


Cne  can  see  that  the  left-haund  side  is  a  vector  of 
the  trial  solution  evaluated  at  tne  collocation  points. 

•  4 

Since  the  term  involving  depends  on  doth  i  and  *,  a 

matrix  will  represent  this  term.  The  b,(t)  will  naturally 

0 

be  a  vector.  If  we  let  a  single  bar  represent  a  vector  and 
a  double  bar  a  matrix,  then  the  set  of  eqns,(9)  through(li 
can  be  'written  as  one  equation  in  this  simplified  notation 

y(ty=Q'b(t)  [  iZ ./ 

■wnere  Q.. 

I  O 

The  first  and  second  derivati'.’es  wicn  respect  to  x 
may  then  be  obtained  from  equation  ill): 
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dx“ 


where  C 


dx“"^  I  ,  _  d^x““^  I 

.  •=  — 7- —  I  and  D,  ^ - = — ^ 

jx.  dX^  jX. 


it  is  now  si.T.pie  to  soive  lor  o\tj  tne  vector  01  c..  [t; 

j 

coefricients ,  since  simple  rules  of  matrix  algebra  apply  in 
this  notation.  Equation  (12)  yields 


b  ( t )  =  Q  ■  y  ( t ) 


Substituting  this  into  eqns.(13)  and (14)  gives 


dy(_ 


dx 


Q'^  y(t)=  A  y(t) 


( 15) 


*1  y  ( '^ )  -  .^”1  q  ..  f  ^  ■ 


(16) 


where  A=  C  and  5=  D  Q~^  or  in  terms  of  the  i-th  collo¬ 
cation  point,  after  some  sim.ple  m.atrix  multiplication, 


dy ( X ,  t )  ■  _ 


K+2 


dx 


-  y  A  :  4  y\  X  , ,  t  ,■ 


1 ) 


42.  '  .  .  -  I  ■''*2 

>  “i.J  ■ 


dx 


J=' 
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One  can  therefore  express  the  spatial  ierivatives  in 
terns  of  the  trial  function  at  the  j+2  collocation  points  . 

Equation  (16)  can  be  directly  applied  to  some  diffu¬ 
sion  process  such  as  that  described  by  ecns.(2)  through  (5). 
Equation  (2)  becomes 


dy ( t ) , 
dt 


f(y(t)} 


(19) 


where  the  vector  dy(t)/dt  represents  the  set  of  time  deri¬ 
vatives  at  the  N-^2  collocation  points,  b'e  note  that  there 
is  no  explicit  spatial  (x)  dependence  in  eqn.(19).  3y  assum¬ 
ing  an  explicit  polynomial  in  x,  eqn.(6),  and  evaluating  it 
at  specific  values,  the  collocation  points,  the  problem  is 
reduced  from  a  partial  differential  equation  with  rwo  inde¬ 
pendent  variables  to  a  set  of  ordinary  differential  equa¬ 
tions  , 


One  can  expand  eqn.(19)  into  its  separate  members  to 
solve  the  problem.  For  the  i-th  collocation  point  we  can 
•write,  according  to  eqn.(19), 


for 


i  »  ^  *  3  I  >  •  •h'*‘2  a 
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One  now  considers  the  conditions  at  the  boundairies 
Xj^=0  and  Since  these  values  are  known  from  the 

original  problem  i.e.  eqns.(3)  and  (4)  one  can  subtitute 
the  boundary  conditions  in  eqn4(20)  to  obtain 


pi^  .y(x-.t))-  f(y(Xi,t)) 

J*' 


for  i  =  2,3i4..,N+1. 


One  should  note  that  g(t)  and  h(t)  are  included  so 
that  the  boundary  conditions  may  be  time  dependent  in  which 
case  g(t)  and  h(t)  would  be  known  or  easily  obtained.  The 
diffusion  problem  has  been  transformed  into  a  set  of  TJ 
simultaneous  first  order  ordinary  differential  equations 
in  N  unknov/ns  with  initial  conditions  and  can  be  solved 
numerically  by  a  number  of  standard  methods  on  a  computer. 

Evaluation  of  the  current  at  an  electrode,  i=  nFAD|^!  , 

ix=0 

can  be  obtained  directly  from  eqn.(17)  once  the  concentra¬ 
tion  profile  is  obtained. 


APPENDIX  D 


COIVFUTS: 


PROGRAMS  USING  ORTHOGONAL  COLLOCATION 


METHC: 


1 
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¥'J 

C 

c 

THIS  PwQijaA.''.  oUL72S  E.i'JATIO.iJ  (50)  TiiSJJjH  (55)  n'li.i  lau 

2 

c 

DJITIAL  JOiiDIXiONo  JSIMG  SPLINE  CGLLOCATIGN  rtSThOG  FOA  T  :i  i 

c 

c 

CASE  3F  THE  ELECTSOL7TS  CQNl’AININJ  NO  CAD.-IIU;!  ION. 

'Z 

^  c 

EQUATIONS  (50)  AND  (51)  AiE  DISCazriESD  FIAST-OAOEn 

c 

c 

OSDINA3I  DIFFEaSNTIAL  EQUAIIONS  «HICa  ASE  SoL/ZO  BY  AN 

c 

SSai-I.IPLiZir  INTEGSATION  S'JBROUTINS  CALLED  DIFSOB  JSING 

c 

1  c 

GEAR'S  aOOriNE.  EQUATIONS  (52)  THAoUGil  (53)  A3Z  SOLVED 

SIJULI ANEOUSL  Y  BY  A  SUBriO'JTINE  CALLED  ESYSia  'NlllCH  CCaZS 

'2 

V. 

F?Oa  lES..  SJjROUriNE  PACKAGE  10  OBTAIN  TiiE  SUaEACE 

C 

J 

CD;!CE;iTaAr..O:ij  of  nitfatz  and  hydboaI  ions.  SQJArioi.s  (54j 

C 

c 

AND  (5  5)  ASi  InZia  BULK  VALJZS.  TEE  CONCENTF  ATION  PriCEI.-ES 

2 

c 

THUS  ODIAINLO  ARE  THEN  USED  TO  CALCU..ATE  THE  CUFREN: 

C 

c 

DSNSIIY  jY  -^UATION  (5b).  THE  CUFaE..:  DENSITY  DATA  AFE 

u 

L 

THEN  USED  TO  FIT  AITH  THE  E 1 P £3 Id ENT AL  DATA  TO  IDENTIFY 

2 

THE  CATHODIC  ANO  ANODIC  aEACIION  SATE  CONSTANTS  OF  THE 

c 

1*^ 

SURFACE  aZACriuN. 

w 

c 

V. 

c 

-  NOaSMCLA.  UBZ - 

c 

s* 

:i.(=;.U;iOER  jF  points  in  X  OlaZCTIOa. 

,,  ^ 

N=NUdBER  OF  INTZaiOa  COLLCCAflON  POINTS. 

c 

n 

NT=NJM3Ea  OF  TOTAL  COLLOCATION  POINTS  INCL’JDINo  TNO 

c 

c 

JCUNDitaY  POxNTS. 

n 

s. 

Y  (1  ,  1)  ,  Y  ( 1  ,  T)  , .  .  .  i  ( 1  ,  ;.•)  =CONCENT3AriONS  OF  NITRArZ  ION  AT 

2 

r 

N  iNTSaiOK  COLLOCATION  POINTS. 

2 

Y  ( 1 ,  :)♦  1 )  ,I  (  1 ,  N42)  ,  .  .  -  Y  (1 , 2*N)  =CONCEN  TR AT  10 NS  OF  H  YOaUai 

z 

c 

UN  AT  N  INTEalOS  COLLOCAIiON 

c 

yj 

POINTS. 

c 

Cl  {  1)  ,  Cl  (2)  -cl  (NX)  =CONCZNTEATICN  OF  NITRATE  ION  AT  EACH 

c 

J 

POINT  IN  X  DISECTIO;;. 

c 

w 

C2  (  1)  ,  w2  (2)  ,  .  .C  2  (NX)  =C0NCENT3ATI0N  OF  H  YDFOXY  ION  AT  tACH 

c 

z 

POINT  IN  X  DiaSCTIUN. 

c 

\  z 

X  3  (  1 )  =  SU  a  ?  .kC  i  CO  NC  E  N  TRATI  ON  OF  NIIXAIE  lON. 

2 

c 

X3  (E  )  =  Su  a  F  AC  .j  CO  N'C  E  N  TE  AT  I  ON  oF  HYOaOXY  ION. 

2 

1  c 

C1J=jULK  CONC ENTFATION  0?  NITaAIS  ION, 

2 

t 

,2B=3ULF.  C  JNC  ENTFATI CN  OF  iJ'YOjOCY  ION. 

c 

DaAX  =  OIFFJSION  THICKNESS,  Cd. 

C 

1 

DS  =  LOCAIIOU  Of  I'HE  SPLINE  POINT,  0  <  03  <  1  . 

c 

i  ^ 

DZS  =  i.NC3EdENT  0?  THE  SPLINE  POINT,  ES. 

2 

ri.iAX=;uxidu;i  .i;iE  re  be  iniegeateo,  sec. 

C 

JT  Id E=  Tldi,  IN  CaZHENT,  SEC. 

vj 

TAJH.iX  =  OIdE.IS  lo'UESS  EAXIdUH  TIdE  To  BE  INTEGaATEO. 

1 

;• 

c 

..  .  .V  J ^  J  1 .1 M  J  1  »« L  J.I S  r  I 2  I  wC  2  2 j'  i 

2 

J 

A a  —  ;v  c,A  2  A  i  u  .i  Ci'  iOa. 

2 

1 

0  4  T  A=  a  i  ACT  i.  j .(  u?J"?  OF  NITEATE  ION. 

(  '-* 

COa^  —  J.J  akiAs^^XO*'!  aAi.w  X>*  A.iOJi.'.- 

J 

Z  ww*  Xu  >«• 

z 

w  0 .1  1  =H  a  r  j  a  0  S  a  .<  wO  US'  AL’ACTICN  NATE  CO.iSTAaT  IN  CATHJD.C 

Ua.v  wCa.  Xw  .*« 

CD^CUaaEa'T  DE.tSTT'Y,  AdP/Cd**a. 

-i  . 
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THE  COEFFICIZ:iT j  -F  Idi  ..’irRATE  A.O  HM.UaI' 

IONS  .\22  1.^JJE-5  \N0  5.26JS-5C0^/S2.J,2E3?ZCTIVE0I. 

«*•«««««*«»> 


C 


c 

«c 


r;iPLicir  oSal*s  (A-h, .^-2) 

OiaZNSiJN  A  (2  0 ,2  J)  ,0  (20  ,2  j)  ,01(00)  ,C2  (60)  ,  C  (20j  ,  K  (3 ,20)  ,  /  /.AA  (20)  , 

1  OIr 1 (20)  ,DIF2  (2  0)  ,D;f J  (2  0)  , ROOT  (20)  , 7ZCT  (20)  , ERRCO  (20)  , 

2  SAV2(2J,20)  ,?^t400)  ,XS(2)  ,AX(1)  ,»A(J)  ,3A5  (1) 

EXTERNAL  A J 

CG.1.10N/L1/CJZF1  ,CjZF5 
Ca220N/L2/20£F2 
CCiA.lON/LJ/N  ,  N  r 
cu;i'io;i/L-*/A ,  D 
COONON/Lb/XS 
C002ON/Lo/0OiiF4 


0  U  tiN  ON/  4-  7  /  .1  S  j 


Cu.i:i  JN/Ld/O  J-F  3,CONl  ,CO  N2 ,0  A2A  ,  JOTA 
cg;i;!cn/lo/y 


INPJ.  DATA  -■ 


RZAJ  (5,3  10)  I'lD,  ;i,  NO,  N1  ,  AL,  3E,  2?S,  TIG  A  X,  IG,  DOS 
REnJ  (5,320)  Cl3,C23,D/iAX,Dr4..'!E,NC,2S 

READ  (5,330)  4,A,-iA,  B  ETA,  CC:i2,GO.i  1 

BEAD  (5, 340)  uSIG,NE2,Z? 

- CALCULATE  SO;i£  CONSTANTS  TO  oZ  USED  IN  T.IZ  F50SBA2 - 

COZ F 1  =  (  1 .  7  02  J- 5)  / 5- 2bD -5 

C0i^42  —  C44S/C  1  5 

CU2F3=  (  1,  J02D-5)  *  Cl  S/D 'U  X/ IS 

CU£F4  =  2,  *  (  1.  902D-5)  *01  3*^6437. /o;iA.{/2S 

COijr5“l4/4jS/  440 

Da=0>''.AX/  <NX-1) 

DTAU  =  (d,  2oJ-5)  *DTIXE/0.'5AX/DNA.v 
TA0i1AX=  (5.  26o-5)  *T  12  AX/DO  AX/D.lAX 
SDTA  a=oTAu 
DETA  =  1./  (NX-1) 

NT=N  •‘..•ia+N  1 
N£;3=2*N 

------  SSI  INITIAL  DI  2  ZNS  I GN 14Z  So  11.12  =  0.  ------- 

T  A  J=  0. 

- DSFIi.Z  7NL  IIAaIN'J.I  A..J  .UNISJ.'.  D  I.'l  ZNS  I  J  N  LSS  o 

4  X  .1  o  I N*  C  .4  2k  N  4  4  0  00  j  o  o  o  1,  N  4 .  i  *1  ^  N  4  o  o  41  .t  4  4  44 .. 

SODooUTiNo  D4.roU5  —  —  ——  —  — 

0  T  .1 0  .'!  X  =  4  0  0  1  *  L»  I A  U 
DrrtJ.1A=2.0*  OTA  J 
«  S I T  Z  ( tj ,  j  j  0  ) 


1 


n  n  r.  r>  n  o  n  n  n  n  r,  n  n  n  n  n  n  n 
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- EVAI.JATE  THE  COLLOCATION  ?0  IN  TS  (BOOTS )  OF  THE 

JACOoI  POLYNOMIAL  CT  OLOEA  N ,  AS  4  EIL  AS  Tr.E 
FlaST  A;0  SECOND  DEBIVATIVES  OF  IHS  POLYNOHIAL 
Hi  TaE  ROOTS.  - 


CALL  JCOtii  (WD,  N,  N'O  ,N  1,  AL,  3S,DIF1  ,DIF2  ,  DIF3,  ROOT) 

WRITE  (6 ,3bOj  (ROOT  (J)  ,  J=  1  ,.iT) 

-  CALCULATE  THE  D ISCR ET IZA TICS  COEFFICIENT  MATRIi  A  - 

10  =  1 

DO  20  1=1, NT 

CALL  OFU?  d  (ND,  N  ,  NC  ,:n  ,I ,  ID,  01?  1 ,  01:2,  0IF3,  ROOT  ,  '/ECI; 

DO  10  J=  1 , NT 
10  A (I, J) =7ECT  (J) 

20  CONTINUE 

- CA-C'JLATS  THE  DISCRETXEATiON  COEFFICIENT  3 - 

ID=2 

DO  40  I  =  1  ,  NT 

CALL  DFOPrt  (;I0,N  ,  NO, Ml  ,I,iO,OI?1 , 0IF2 , 01:3  ,  ROOT  ,  V  £C„  ) 

DO  30  J  =  1 ,  NT 
30  3  (1,  JJ  =VECT  (J) 

40  CONTINUE 

-  INPUT  THE  VALUES  OF  ARGUMENTS  TO  5£  USED  IN  THE 

INTEGRATION  SU'dSOUTINS  DIFSU3 - 

HF=1 

JSTARI=0 
MAXDEE  =  7 
DO  4  5  1=1, N  EB 
THAI(I)  =1.  J 
45  CONTINUE 
1  =  0 

-  INPUT  THE  INITIAL  CONDITIONS - 

DO  5  0  J= 1 ,  N 

Y  (1  ,0)  =1. 

Y  (  1  ,  J  ♦■N)  =CO&  F2 
50  CONTINUE 

XS(  1)  =  1- 

\  mJ  ^2  }  . 

W  R II  E  (  o ,  3  s  0  ) 

bO  riHE=IA  J  =  D.i  Aa  »UMA\/  (5.  26D-5) 

nfciD  PRlN«  lHci  CJRRciNi.  DitNoITi  AT  i.HIb 
CCNCENTSATi  JNS  OF  NITRATE  A..0  .iiLROIi  lo  ..S 
COLLOCATION  POINTS  — — 

CALL  CJRENT  iCD, A, Y) 


*  Ciii^^iiLA^i 

PalN.  THE 
1%  L  .  li  E  N  T 
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c 

c 

c 


t- 


c 

c 

c 

c 


TlMEfCDflCJ(t)  /  .C^(2)  /  KALACj^IXZ^kifM^ 
iiair£(3,7JJ)  ('/  (1,L)  ,L=l/.i2A) 

-  iVrtuJATE  A:JD  ?iii:JT  r;i2  COiSCZl-JTP.AnC'iS  Or 

:iir.-<AT£  AMC  HYOaOXY  IQM5  AT  2ACM  POIUI 
IN  1  OIRSCTIOU  AT  i'liIS  11.12 - 

CALL  CONCEU  (C  1,C2,  NX  ,N  D  ,0  ET -A,  ROOT,  DI?  1  ,'i) 

.iaiTE  (6, 230;  (C 1  (  J)  , 0  =  1 ,  N  X) 

WRITE  (6, 2-jO;  (C2  (  J)  ,J  =  1  ,  NX) 

- CE2CX  IF  THE  I N  TEC  A  Ai  I  ON  TI.IE  AT  THIS  .lO.lENT  'iE\CnZ6 

I  ii  L  ilAXHIUH  TI.NE  TO  32  INI2GBATED  — - - 

IF  (TAU.  OE.  T  AO  .1AX)  GO  TO  ^  jO 
55  !=!♦  1 

IF  (I.  EO.  1 )  Gu  TO  9  0 


V. 

C 

c 


TEST  IF  THE  CO NCE NT aA TI ON  C?  NITSATE  ION  AT  THE  NIH 
COLLOCATION  POINT, IE,  T.42  LAST  ONE  3EF0SE  THE  SPLINE 
POI.iT  IS  WITHIN  O.CGOl  OxHENS  lO.i  LESS  CONCENTSATION 
JNir  OF  THE  30UNDAAY  CONOiriJ.i. - 


C 

c 

c 


IF  {  (  1. -Y  I  1,  j;  )  .  LT.  0.  0001)  NO  TO  30 

-  IF  THE  TEST  FAILS,  THE  SPLINE  POINT  *3  INCI.iASEJ, 

;iOVEO  FUPTHE?  INTO  THE  SOLJTION. - 

CALL  CHANGE  (SS,  OZS,  COEF3,CijEF4,COZF5,  FA  CTO  H,  0TA.J,0rAJ.'1I, 
1  J  START) 


0  I A  0  .'1 A  , 


EVALUATE  THE  CO  NCEN  laATIONS  OF  NIiaAl'i  AND  HYOAOxY 
X  J  .4  S  AT  T  .1  c  COLLOCATI'jN  POINTS  AX  THx  NaW  CCO  a  D  Z.i  a  4.' z 
OJE  TO  THE  INCREASE  0?  THE  SPLINE  POINT.  - 


C 

CALL  EaPANO  lY,SOOT, 0IF1,FACT0  5, NO) 

C 

C  -  liiiiOSATE  EQUATIONS  (50;  AND  ;51)  OSINj  THE  CONCENTFATIu 

C  li  4  Tnr.  COLi.OCATxON  POxN.S  AT  xHIa  TI.IX  TO  O&lAx.*  T.IE 

C  COnCENTP  AT  IONS  AT  THE  SJCC  ESS  IVE  11 'I  - 


9u 


C.iLu  0  IF S  0 i/  ( N  EH  ,  r A  0 ,  Y,  SAVE,  OTA  J  ,  OT  AON  I  ,  0 T A  0  HA 

X  a  a  0  5 ,  K  F  L  A  G ,  0  s  I A  a  r ,  .1 A  X  0  5  a ,  ?  ’ ) 

IF  (K7LAS,  LT.  1)  GO  TO  190 
jO  1  *LiTA0 
JTAU.U  =  2.  0*0 TAG 
■.F{A.Gx.xj  jU  to  9  3 

IF  U  (1  ,  1 )  .0  1.  1 ..  OP  .  Y  (1,2)  .or,  1 - ..R  .  Y  {  1,3;  .  GI. 


P  X  ,  .1  r  ,  .  .1  A  , 


•  0,^.t^t,x;  •  ji.i. 


no  on  c  i 


1 


lii.2 


.  JE.  DAoi  (  i  ( 1  *  5)  -  1_  )  .  U  r.  0.  J304)  ju  ro  '90 
GO  TO  95 
9  3  ;n=i/ic: 

1*10 

IF  (I.  E  .12.  OE.  TA:J.  GE.TA'JMAa)  GO  TO  95 
GO  TO  tJb 

- CALCULATE  THE  SURFACE  COJ«C EG T  HATIO NS - 

95  XX(1)=Y(1,1J 
ITSaN=5U0 

CALL  ISiST.l  (AUX2,Ht,NSIG,  NEZ ,  XX  ,1  i’ ESN  ,  «  A  ,  B  A  E,  lER) 

XS  (1)  =  .CX  {  1  ) 

XS  (2)  =-i  -  »CJEF1*  (A  (  1 ,  1)  '"XX  (1)  *-A(  1,  NT)  «1. 0)  /A  (  1  ,  1)  -  A  (  1  ,  .'il)  /A  (1  ,1 J 
1  *CJt?2 

DO  150  J=1^G 

XS  (2)  =  XS  (2)  -  3.  •CCEF1  *A(  1  ,  J*- 1)  ♦X  I  I ,  J;  /  A  (1  ,  1)  -  A  (  1 ,  J+1) 

1  •X  (  1,  J  +  N)/A  {1,  1) 

150  co.'rriNUL 

GO  TO  60 

190  waiIS(tj,  240)  DX,SDTAU,  NX,DTAU,DTi;}E,01AX,  AL,3  5,  ZS  ,  EPS,  XFLA  J, 

1  e?,ns:g 

XnlTE  (o,  2>*5)  BETA,  GAMA  ,  CO  ;11  ,CJN2 

-  FDRIATS  FOR  INPUT  AGO  OUTPUT  ST.ATE1SNTS - 

2  30  FOR3AI  (/,25X,  ’NOS-  1X,9D  1  1.4/30X,  9D1  1.4/30X,  9011. 4/30X,  90  n.  4 

1  /30X ,901 1. 4/30X, 901 1.4/30i, 901 1. 4/30X, ^01 1.4/30X, ^011.4 

2  /J OX, 901 1. 4/3 OX, 90 1 1-4/3 OX,  9011. 4/3  OX, jO 1  1.4/30  a, 90  11.4) 

240  FOR1AT  (//,  bX,  '  X  I  N  T£  R  V  A  L=  '  ,  3 1  9 , 3  2  X ,  '  I  NITI  AL  TAU  STEP  SiZ£=' 

1  , 01 2. 4/5X , ' NO  OF  POINTS  IN  X  DI?=  ',I3,29X, 

2  'LrtST  TAU  STEP  SIGE=  '  ,  0 1 2 . 4/5  X ,  '  RE  AL  THE  IUTSayAL=  ' 

3  ,  J10.4,24.X,  'DIFUSIOX  LENGrH=  '  ,  0  1  0- 4//5  X ,  '  AL=  ',:-7.2,30X, 

4  •uE=  • ,?7. 2,Z5X,'iS=  • , ?7 .4//5X, • E?5=  ',D10.2,27X, 

5  'r;FL.»G=  '  ,  I3//5X,  '  E?  =  •  ,  0 1  0 . 2 , 27  X  ,  ’  N  5 1  G=  ',lj) 

245  FUtilAI  (//,5X,  'REACTION  ORDER  OF  .403-  ION  =',F7.2,25X, 

1  ’aSACriON  ORDER  OF  OH-  ION  = ' , F7. 2//5 X ,  '  1  ST  CONSTANT  =' 

2  ,  j12.4,25;(,  'END  CONSTANT  =',012-4) 

2  50  FOol AT  I// 25  A, 'OH- ',21,9011. 4/3  OX, 9011. 4/3CX,90  1  1 .4/30X, aO 11.4 

1  /JOa,  :»0  1  1  .  4/3  OX,  90  1  1.4/3  OX,  90  1  1. 4/301,  9  0  11.  4/ JO  X  ,  9  D  1  1 . 

2  /30a,901 1. 4/30X, 901 1.4/3 OX, 90 1 1 . 4/3 0 X, jO 1 1.4/30X, 1 1 .4) 
310  F05.1AT(4I3, 2F5.  1 , 2  0 1  3.  2  ,I  5  ,  Fd,  4 ) 

3  20  FOa.'lAT  (4012.5,17,  F7. 4) 

3  30  FOS;i  aT  (29  7.  X,  20  1  2- 4) 

3  4  9  FOa.1  AT  (2X5,  01  2.  4) 

3  50  FO  3.'!  A i  1 1  5  X  ,  1  0  F  1  0  .  6  ) 

3  60  FOa.l  AT  (//,  01-.4,1X,D12.4,2.v,2DjJ.1J,-15,F2  0.7j 
J7  0  FuoHAl  (/,  2a,oF2  0.  6/2  X,  oFJJ.  o) 

3  90  FQ  R.1  AT  (2  X, 'COLLOCATION  POI..  IS;  •) 

3  90  FOR.l  at  I//,  3  A,  'T:  A.S  '  ,  12  X  ,'CJ3R£Nr'  ,45X,  '  ION  CO.iCEMTL.MTv,.'.  '  ) 

7  00  FOR, -I  AT  (3  OX,  10F1  J.  5/3  Ji  ,  1071  3,  5) 

STOP 

END 
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SU55  0UriN  Z  JCU3I  (N  C,  M,  NO,  ;j  I  3Z,  JI?  1  ,  Di:-  2,  JIF  3  , 1'.OCl) 
I.'lPLICII  F  L  AL*  d  (  A  -  H,  0- 0  ) 

OI.IiNSiON  DiF  1  (:iD)  ,DIF2  {ND}  ,OIF  J  ^NO)  ,  FOOT  (ND) 

S03R0UTINt  cVALUATSo  THE  COLLCCATICN  POINTS  iSOOTo)  OF  THE 
JACOBI  POLINOJIAZ  OF  05D£B  N  AS  ’nELL  Ai  THE  FIRST,  SECoNO 
AND  THiaO  derivatives  OF  THE  POi-YNOJIAL  AT  THE  ROOTS. 

- IN  PUT  PAR AJETERS - 

INTEOEr  ND=THE  DIJENSICN  of  vectors  DIFI , DIF2,DI?3 , ROOT. 
INTEoER  N  =  ihE  DESS  EE  0?  THE  JACOBI  POIV  NOilI  AL ,  I.E.,  THE 
NOSSER  OF  INTERIOR  COLLOCATION  POINTS. 

INTEGER  NO=D£CIDES  WHETHER  X=0  IS  INCLdDED  AS  A  COLLOCATION 
POIrT.  no  HOST  3E  SET  E^dAL  TO  1  (INCLJDING  X=Oj 
OR  0  (EXCLUDING  THIS  POINT). 

INTEGER  Nl=hS  FOR  NO,  BdT  FOR  THE  POINT  1=1. 

REAL  A..,3S=IIiE  PARAMETERS  OF  THE  d.iC03I  POLYNOMIAL. 

-  OdTPdT  PARAJETERS  - 

REAL  ARRAY  EDO r=ON E- DI ME N SI  ON AL  VECTOR  CONTalMING  ON  EaI: 

THE  N*N0«-N1  EEROS  OF  THE  POLYNOMIAL  USED 
IN  THE  COLLOCATION  ROJTINS. 

REAL  ARRAY  =rKR2E  ON E-DIMENSIONAL  VECTORS  CONTAINING 

DIF1,DIf2, DIF j  ON  EXIT  THE  FIRST,  SECOND,  AND  THIRD 

DERIVATIVES  OF  THE  POLYNOMIAL  AT  THE  EEROS. 


C  - FIRST  ..VALUATION  0?  COEFFICIENTS  IN  RECdHSION  FJtvadLAS  - 

C  RECURSION  COEFFICIENTS  ARE  ST03SD  IN  DIF1  AND  DIF2  ~ 

C 

A  9=  A  B  £ 

AD=BE-AL 
A?=3E*AL 

DiFI  l1)  =  iAD/(AB*’2)  ♦•1)/2 
DIF2  (T) =0 

IF  (N  .LT.  2)  SO  TO  15 
DO  12  1=2, N 
^1=1-1 
-=A3  ♦•2''-1 

DiFI  (I)  =  (AS  i>A  D/2/  (2*2)  ♦  1)  /2 
IF  (I .  NE.  2)  GO  TO  1  D 
DIF2  (I)  =  iAj*rt?*’E1)  /Z/Z /  (E*’1  J 
S u  TO  12 
2  =  E*  E 

Y=21  *  (A3»u1  ) 


10 


c 

c 


'i=i* 

0x72  (I) 

12  COlil  IN  J  E 

— - -  3J^i'  0£7  27  .'I  I  NA  TI  0  N  32  NZmI^N  .'12iT0i>  'nir.’i 

OF  P2i\/I0aSLY  DET2a.-!I.,£;j  auOTG  - 

15  X=0 

Du  3  5  1=1  ,  N 
20  .<D=0 
XN=1 

;(Di=o 

.<  N  1  =  0 

uO  2  X  J  =  1  f  ii 

X?=  (Oir’1  (J;  -x)  ^XN-CIF2  (J)  *2^ 

:<?  1=  (DIF  1  (J)  -  .<)  *  aN  1-Dir2  (  J)  *2D1  -  CN 
XD=aN 
XD  1=XN  1 
XN  =  XP 

22  \;n=x?i 

2C=1 

2=x;;/>;n  i 

IF  (1 .  2*.  1 )  DU  :0  3  3 
DO  25  J=2,I 

25  <.C  =  2C-2/ (.(-FOOT  (J- 1)  ) 

30  2=2/20 
X=  X-  - 

IF(DABS(-)  .  OX,  1.0-09)00  TO  20 
sour  (I)  =X 
X  =  X*.  00  3  1 
35  CUNIINOZ 

- ADj  OVFNIUAL  CGLECCA'XION  POixrS  AI  X=J  03  A=1 - 

NT  =  N  +N0  ♦•N  1 

\jO  ro  *^5 

Do  40  1=1,  j 

J  =  .|4  l-I 

40  aour  (J 4 1 )  =20ux  (j) 
nOOI  (1)  =0 

45  If  (N  1.  2^.  1  )  JOOr  (NT  )=  1 

NJa  2VALJATE  DEFIVATI/ZJ  UF  ?UL  I N03I AL - 

Do  5  0  I  =  1  ,  ;i  ? 

>:  =  r.JoI  (I; 

DIrl  (I, =  1 
DIx-2  ll;  =  J 
JlF3  (I) =0 
DO  0  J  J  =  1  ,  ,1  i 

,0  3  0 

'i-\~  .lOo  r  ( J) 

01 F3  (1)  =F*JIi-  J  (I)  4  3*0172(1; 

01:2  (1/  ~L*jLt2  (I)  42*jirl  (Z) 


1 


fj  n  (j 


14 


3IF1  (!)  =1*3  ii’  1(1) 
50  OCNli.NJi 
FZT'J  a:i 
e::d 


SUoaOUriwZ  ur  ops  (N  D,  :10,  N  1  ,  I,  13,  Oif  1  ,  3IF2, 3IP3  ,  SOOT,  ViCT) 
I.1PLICI7  r<ZAL*d  (A-H,0- 3) 

Dl:iENSION  3IF1  (S3)  ,31x2  (HD)  ,DIr  3  (:i  3)  ,auOT  (UD)  ,  7  EC  r  (SO; 


C  cajciouriss  EVA^OATZS  DISCSZTICAriCS  coeffiliest  .iaiaicec  c 

<:  AS3  GAU53i.AS  vUAOSATURE  WEIOilTC,  SG.1ALICE3  TO  S  J  .1  1-  C 

c  c 

C  - IS  pur  PAAA.'IETEES -  C 

c  c 

C  ISrECER  S,  JO,  S1=AS  IS  S'JBSOUTIME  CCCOI.  C 

C  ISTECES  I  =  rHi  INDEX  CF  THE  NODS  FOj  JHICH  THE  JEIGHiJ  A3Z  v. 

L'O  aE  CALCULATED-  C 

INTEGER  1U=1N 3ICAT0S-  13=1  GIVE3  THE  JEIGHT3  FO 3  3T/3X,  C 

13=2  GIVES  THE  WEIGHTS  -03  321/3X2,  AND  13=3  w 

GIVES  THE  GAUSSIAN  WEIGHTS.  THE  VALUE  OF  I  IS  C 

C  Ia:vx,LEVANT  IS  THE  LAST  CASE.  C 

c  r.  E  A  L  A  3  ?  A I  a  0  J  r ,  c 

c  JlF1,3IP2,3iF  3  =T;IE  ONE- 01  HE  NSIOSAL  VECTORS  COMPUTED  iJ  C 

C  3U3KOUTISE  JC03I.  C 

c  c 

C  - OJTPUi  PASA.ISTSSS -  C 


C  PEAL  AflBAI  VECT=THE  COMPUTED  VECT03  OF  -EIGHTS.  C 

c  c 

NT=N  *N0  *'3  1 

IF  (I  3,  E;.  j)  Gu  TO  1  0 
30  20  J=1,Lii. 

Ix'(J.JE.I)GU  TO  21 
IF  (I  3.  .nE.  1 )  Gu  TO  5 
VECT  (I)  =3  I?  2  (  I)  /Cl  FI  (I)  /2 
Gu  TO  20 

5  VECT  (I)  =3xF  J{  I) /DI  FI  (I) /3 
G  u  . u  20 


21 

T=  aO  uT  ( I ) 

-auUT  (C) 

VEC:  (J)  =3 

IF  1  (I)  /  3  1 1-1  (J)  /T 

xr  |I3.X|^. 

2)  vEcr(j)  =vEcr  (C)  *  {3IF2  (D  /jifi  ^i 

r  > 

CuMiSUE 

u  3  .0  3  3 

10 

1  =  0 

Ou  2  5  J  =  1 

,G  r 

::=53uT  ij) 

A  X  =  X  •  I  1  -  A ) 

I F  (  N  0 .  E  «  .  3  I  A  *  =  \  X  /  X  /  X 

.X  (ii  l..^^«0)  iiX  =  Aa/  (1*X)/(1“a) 


1-^6 


VZCi  (J)  =AX/ JIF  1  (J)  **2 
25  f  =  Y  +  /  E  C  r  i  o  ) 

CO  0  0  J  =  1 ,  X  1 
60  VECI  (J)  =  VnCT  (J) /Y 
50  HETUaN 
END 


s  jaaoariNE  i.itb?  (n  c,  nt,  i,  aooT,  2i?i ,  xi:itp) 
i/i?LlEIT  EEaE*3  (A-H,0-2) 

d1i.en51un  douT  (nd)  ,difi  (nd)  ,;{iNr?  (nd) 


c 

C  SOOaOUTINE  EVALUAIES  THE  LAGN  A.'.’Ol  AN  I  NT 


C 

EDP'jLAriOS  OuEFf  lOIZNTi 


C  - INPUT  ?Aii  ANZTZaS -  C 

c  L 

C  INTZOZa  Nr  =  rnZ  total  N'JlijED  OF  COLLOCATION  POINTS  (  =  N  ♦■:.0  + ..  i )  C 

C  FOfl  WHICH  THE  VALUE  OF  I'nZ  DEPENDENT  VASIALLE  Y  C 

C  IE  KNOWN,  C 

C  HEAL  A-TtiE  A3ECISSA  X  WHZZE  Y  {X)  IJ  DESILEO.  C 

v-PiALABuAY  C 

C  ROOT,DI?1  =ZDLLOCATICN  POINTS  AND  DESIVnTIVES  OF  NODE  C 

C  PJ^YNOHIAL,  DEaiVED  IN  jUSaoUTINZ  JCo3I.  C 

c  c 

C  - OUTPUT  PAtiAHSTEaS -  C 

•I  C 

C  SEAL  ASBAY  XINTP*THS  VECTOB  OF  INTSaPOLAIION  WEIOKTS.  C 

>1  C 


'♦♦***  ♦  I 


:  «  «  i|C  «  « 


?CL=  1 

DO  5  I  =  1  ,  N  T 
Y=X-  BOOT  iD 
\INTP  U)  =0 

IF  ^Y.  Es^.  0,  DO)  XINTP  (I)  =  1 
5  POL=POL*Y 

IF  l?OL.  E,;.  0- DO)  GO  TO  10 
DO  o  1=1, NT 

o  NlNT.’(i)  =POL/Diri  (I)/(I-3Uo:  tl)  ) 
10  AETURN 
END 


G  U  J  n  o  U  T  ^  4i 

1 

ToP-icir 

Jl.-iZNSIDN 

3 

D  ,  C  all  4*  H  i.  J  iN 


E  Dxa'ajUo  (4^,  L,  ,,'JA  /  4J,i1,>l>.X4i,.‘4  4'liAX,4rfP>J,  .'t  £  ,  £ 

•JSTA  R  r,  UaAXDER  ,  P  .1) 

(A  *44,  ^”*4) 

Y  to, 20)  ,'Y.1AX(2J)  ,  SAVZ(2'*,tO)  ,EI\liOat2J) 
A  to)  ,PERT3T(7,2,j) 

X  txO)  ,?JL{20  ,20)  ,,Ct20) 


.1a 


t 


ki  <m/  I*' 


,2ti  k-t  0  J  ;  , 


^  h  Li  fA  VO  ^ 


o  n. 
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\ 


i 


I 


c 


3  i  i  d  S  'J  J  C 1  It  3  X  in  X  li  tj  ?  A  T 
CSDaS  i;^UAriON5  A3 
0 Z  31  EP  OF 
DSC3  3A3E0  .liiaX.'i  TH 
AS  EaaOE  A  STEi;'  AS 
STEP  SHiiOfi  HU  iC 'd  IS 
EACH  CUilPJijENr  07  7 


S3  A  SZT  OF  .<  OTiOINASl  DiFFSAE.'iXIAa  FI33X 
0ESC3IJSD  It!  i^JAriOliS  i^O)  A.iO  (51)  OVxc 
a  AT  SACa  CAaL,  a  HAX  HE  XNCREASEO  OF 
S  AAilSS  :1.UN  TO  a,iA.C  IS  OSDSA  TO  A'-Six/S 
POSSIBLE  *:iIL£  itUT  CuaiUlTINu  A  SINGLE 
LA  no  Z  3  THAN  EPS  IS  TSS  L~C  NO?.''.,  NSaLE 
HE  Z330?  IS  DIVIDZJ  Bi  THE  CO  .1?  Ji.E  i.TS  uF 


Ya  Aa- 


THS  PaoGHA.i  calls  FOL'H  SJ2RUUTISES  .SA.1EJ 
OIFfS.'J  (I,  Y  ,0Y) 

HATiNV  (P  J) 

SOLVE  ( it  ii ,  p 'J  L ,  B ,  K) 

PEDES  V  (7,  Y,?;i,.a) 

THE  rIHST,DIrFUa  , EVALUATES  THE  DERIVATIVES  OF  ME  DEPEinxi.: 
VASIAdLZS  STOSED  I.V  Y(1,I)  70S  1=1  XU  :i,  AND  3T03ZS  T  .lE 

DESIVATIVES  IS  THE  ASSAY  CY-  THE  SECOaO  IS  CALLED  ONLY  IF 
a  7  IS  SET  TO  1  05  2  703  STI77  .lEI.ioDS.  IT  PESFOS  IS  A  rlaJT 
STAGE  OF  GAUSSIAil  SL  I.a  I M  ATI  Cii  SOUTINE  OF  THE  iJ  3Y  N  E.AXSIa 
STOSEO  IN  THE  A33AY  ?W{M,a).  SOLVE  PSaFOS-IS  THE  aACK 

sussTiTUiioN  Process  of  the  gaussian  elieinaiion  sojtinx. 

pMa^.tV  xo  cjSc.0  ONx.Y  xi*  .*1r  Io  1,  A.IO  xiix  P/^aix^t*.. 

DEPIVATIVES  Of  T.4E  DIFFESENTIAL  i^JATIONi  PSOVI.OEL  BY  OIFFJN. 

- PApA.aETESJ  — - 


N=THE  iiUHBES  OF  FI3ST  GSDcS  GIF? S3 nN  TIAL  EQUATIONS. 
T=ril5  INDEPENDENT  VARIABLE, 


c 

Y= 

A  3 

3Y 

N  ASSAY  CONTAIN 

ING  THE  DEPENDENT  V 

A-HIASLES  And 

C 

:.4e; 

h  33 

ALES  DERIVATIVE 

S.  THE  DEPENDENT  VAEIA3LE3  AhE 

c 

c 

3  A  ^ 

ED  I 

N  1(1,1),  1=1  1 

C  N, 

i'*  < 

VE=A 

BLO 

CK  OF  AT  LEAST 

12 *N  FLOATING  POINT 

c 

c 

i? 

Y  ra 

£.  oUSRO'JTINE. 

c 

a= 

THE 

3  T  3  ? 

SIC-  TO  3£  ATT 

SaPTEO  ON  THE  NEET 

STEP. 

w 

H.a 

IN=r 

HE  iUNiau;'.  STEP  SIE 

E  THAT  '.ILL  BE  USED 

F03  THE 

V- 

I 

*1 T  Z  3 

.NATION, 

c 

a;-! 

A  A—  T 

HE  a 

A  A  X  .a  U  .1  SIZE  TO 

sale  a  THE  STEP  w'lLL 

c 

c 

EPS  =  Til 

Z  *Z3r  CCNSTA 

NT- 

c 

EF 

=TaE 

aUD  IN  01  CATO  3. 

0 

AN  ADA.iS  PaEOICTOS  COaSECTOS  IS 

USED. 

c 

1 

A  aULTI-STEP 

aSTHOD  SUITABLE  70r 

STIFF  S  Y  S  T  E .  i  S 

4. 

J 

IS  USED, 

c 

z 

THE  S  A  a  E  AS  C 

AS  5  1,  S.YC  EPT  THAT 

T  .i  I  S  S  U  3  5  'r  U  T  i  N  - 

C0.1PUTES  THE 

PA  3T lAL  DERIVATIVES 

J  k  xi  ^  Ca  X  w  w 

DIFFER ENCI NG 

OF  THE  OESrVAIIVES- 

c 

r 

'/  •/ 

A.v;  =  A 

.i  A  .a  J  A  i  u  r  V  L  CC  AT  I 

ON  S  li'dLC.i  CO.'iT.AI.«S 

X.  .1  w  ■'  .i  Z  xi  J  .'‘i  Ok 

z 

5 

AC  H 

Y  SEEN  SO  FA?- 

'w 

c 

ETv 

.305= 

AN  A 

BN  A  Y  0  7  N  ELEilE 

3X3  V* 0 T  :i  j.  .i 3 

c 

U  4*  Z 

JCEp  ESSO?  IN  £ 

AC  a  co.i  po.n  E.N 

c 

V  ^ 

Arf  .i  VJ* 

A  CO 

a  P  -  E  r  1 0  CODE- 

z 

'* 

J  3 

T  A  3  r 

=  AN 

I R  ?  R  T  r:  D I  CA  TO  3 

z 

.<  0  E3 

=  rHE 

jE?:v\ 

riY—  r>’iAr  z.ijjzj  31. 

V  a  4.  J  • 

S, 

c 

?-yl 

-A  3 

wU  w  »\ 

OF  rtT  LEAST  N  *  *2  fLO.iTING  pix  I N  T  LOCATIONS. 

z 

1 


1^8 


DAT^  PEHiJi  /2«J#a,5^7.33j,10*4*I,jJ»7,17,15#1.J/ 

1  2.0,  1^.D,24.0,37.b9,33.3J,70.0  3,o7.‘j»7, 

2  3-  0,  Q.  0,  9,  1  o7, 12-5,  1  5.  9d  ,  1-  0, 1 . 0, 

3  12.  0,^4. 0, 3  7. d9,  53-33,  70. Od, 3 7. 9 7, 1 .0, 

4  1.,1.,0.5,0.1667,0.04133,O.OOd2G7,1.0, 

5  1.0, 1.0, 2. 0,1.0, .3157, -07407, .013 9/ 

DATA  A  (2)  /-  1.0/ 

IHET=1 
KfLAO  =  1 

I?  I OSIASI.  LE. 0)  OC  TO  140 


1  00 

DO 

1  10 

1=1,  ;i 

DO  1 

10  0=  1  ,  K 

1  10 

3A  VE ( J  ,  I  j  = 

Y 

(O.I) 

:iULO  = 

ii  3  E  « 

IF 

(d. 

EQ. HOLD) 

GO 

TO  1  30 

1  20 

A  C  'J  li 

=  11/ 3  OLD 

13 

ET  1 

=  1 

GO 

ro 

750 

1  30 

:igoL  D 

=  N  >«! 

TOi-D  =  I 

k\  A  C.  U  ~  1.0 

IF  (uSTAdr.OI.O)  GC  TO  250 
GO  ro  170 

140  IF  ( JSTA.i:- £^. -1 )  GC  TO  loO 
.iQ  =1 

:J3  =  N 

:n  =  .S'*  10 

N2  =  N  1  +  1 

;i4  =  N*»2 

:i5  =  Ni  ♦ 

;i  b  =  li  5  ♦  1 

CALL  DIFFU^I  (T,Y,3A7£  (N2,1)  ) 
00  150  I  =  1,3 

150  1(2.1)  =  3A7E  (N  1 +1 ,  1)  *.l 

;i:(2A  =  i 
K  =  2 
GO  TO  100 

160  IF  (  Nv- E'2»  ll  vO  L  0)  JSTABT  =  1 
T  =  TOLD 
NC  =  N>^OLJ 


K  - 

’  :<i» 

♦  1 

GO 

TO 

120 

1  70 

IF 

(IF 

•  c  0  * 

a  0 

TO 

1  dO 

1  F 

{.ig 

•  0  r  • 

0  ) 

vj  G 

TO 

1  90 

GO 

ro 

(22  1 

■) 

f  *• 

22, 

223 

,224, 

22  5,  22 o)  ^ 

1  30 

LF 

( 

•  J  A  • 

7) 

GO 

TO 

1  90 

ro 

(.211 

12, 

2  13 

,214, 

2  1  5 , 2 1  6 , 2  1  7 )  ,  3 

1  90 

c\r  L 

A  cj 

—  “  0 

F  FZ 

•J3;i 

2  1  1 

All)  = 

- 1. 

•J 

GO 

:u 

2  30 

;  12 

A  t  ■ 

)  “ 

5  J 

JOO 

oooco 
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A  i3)^ 

=  -  3 .  5  uo  0  OOO  0  0 

sj  O  4 

u  3  3  J 

1i 

A  (1) 

=  -J.  -*  1o6’jo6o66ou66h7 

A(J) 

=  -0. 750030000 

A  ('0 

•=  -J.  1  o  o6  6  (1 Q  6  fcb  66  b  67 

>j0  L 

0  330 

2  1U 

A  (1) 

=  -0.375000000 

A(J) 

=  -  0.  9  1  o6 b c6  6  6666  6  66 7 

A  (4) 

=  -0.  3  3333333  33333333 

A  (5) 

=  -0.04l6666666666666o67 

■30  10  2  30 

2  15 

A(1, 

=  -0. j4oo1  11111111111 

A  IJ) 

=  -  1 .  0  4  1  66  6  0  6  66  66  7 

A  (-.) 

=  -0. 4  Jo  1  1  1  1  1  11  1 1  1  111 

A  ij) 

=  -u.  10  4lbo6o6b6666b7 

A  (0) 

=  -0-00o  3333  33333333J33 

'j  0  r 

0  3  3 

1  6 

A(Ti 

=  -0,33'aJb  11111111111 

A  t3) 

=  -1.1  *16066660066 66 7 

A  (4) 

=  -J.o2dOOOOOC 

A  i5) 

=  -u.  177  )6  33  3  3333333  53 

A  (0) 

=  - J. JziOOOOOCO 

Av7)_ 

=  0.  00 IjSOSSSSddBOOiiO 

Go  r 

0  230 

■) 

«■ 

17 

All) 

=  -U.  3 1359 193  12169312 

A  (3) 

=  -  1.  2250000QCO 

A  {*) 

=  -0.  7  5  1351  35  13513519 

A  (5) 

=  -0.  26320o33  3333  33i33 

A  (0) 

=  -0. 04o6 1111111111111 

A  (7) 

=  -  0.  0  04  36  n  1  n  11 1 1  m  1 1 1 

A  (0) 

=  -0.  0  00  19341  26934  12o9d'4 

GO  r  J  23  0 

2 

2  1 

A  (1) 

=  -  1 . 0 J4OOOOO0O 

0  230 

222 

A  (1) 

=  -0.ooo66tj6666oo6o6  7 

A  13) 

=  -u.  3j  03  3333  33  33333  3 

r» 

0  23  ) 

■) 

t. 

23 

A(  1) 

=  -vj.  5454545454  5'*  5-*54  3  45  3 

A  (3) 

=  All) 

A  (4) 

=  -0.o9090909G9090  90;»1 

30  r 

0  230 

T 

24 

A(1) 

=  -0.4o0000Q0C 

A  ( J ) 

=  -j./OoOOOOOO 

A  14) 

=  -O«2o0000000 

1  P) 

=  -U.J20J009JC 

0  03  J 

■) 

<i5 

All) 

=  “  J  .  4  J7  9  5o 2  0  43  79  5  62 

A  (3) 

=  -  ).  ..  2  11^7  43  32  n  o73  J 

A  1  4  ) 

=  -J.  3  luZ  1  6a  7  31  02  1  Jbj 

A  1  5) 

=  -  0 .  J  5-*  7  44  5 2  55 -•  7  4 45  2  c 

— 

-S') 

^  ..--J.  :( 6  3  50  36  4  96  J  3  0  4 

JO 

10  J  Ja 

1 

2  0 

All) 

=  -J.  *vJo1  6  32  6  53  Jc  1225 

A  (J) 

=  •u.4*;0o3  4920o  3492Jo 
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A{4)  =  -0  . -I  1  0  6  6 1)6  5  66  oo  o  o*D  o  7 

A  (5)  =  -0- 0  7^20476  1904761 

A(b)  =  -0.  J  1  1  90476  1904  761  5 
Al7)  =  -  0 .  J  00  566  39  34  240  Jo  2  3  2 
230  X  =  MQ+1 
I  Oolj  B  -  K 
:iTY?  =  (4  -  ;i7)/2 
EMQ2  =  .  5/i-' LOAT  (N'O  ♦  1) 

ENOi  =  -  5/ELJAr  (NO  ♦  2) 

ZNwl  =  .  5/:- LOAT  (NQ) 

PSPo'n  =  E?5 

ZUP  =  (JZaYoi  (NQ,  MZY?,  2)  •?£PSd)  **2 
3  =  (PSaiST  (NO  ,«TY?,  1)  *?£?Sil)  **2 
■£0-<N  =  (PEoXiT  (NO,  HTY?,  3)  •PEP5H)  »*2 
IP  t2D«M.£2*0)  GO  TO  730 
BNO  =  £PE*iN03/DFLCAr(N) 

2  40  INEVAL  =  ;ir 

GO  ro  (250,od0),IaEr 

2  50  r  =  r  «■  ti 

DO  2b0  J  =  2,X 
DO  260  J1  =  J,K 

J2=X-J1  ♦  J-  1 
DO  2cJ  I  =  1,N 

260  1  (02,1)  =  Y(J2,I)  ♦  Y(J2*1,I) 

DU  2  70  I  =  1,  N 
270  Ea£i3a(L)  =  0.0 
DO  430  L  =  1,3 

CALL  DIF PUN  ( T , Y ,S A  VS  (  N2 , 1 ) ) 

IF  ( I*  SVAL.  LT.  1)  GO  TO  350 
IF  (.'IF.  £2. 2)  GO  TO  310 
CALL  PEDE3V  (T,Y,  F4,H3) 

A  =  A  (  1 )  *  3 
DC  230  I  =  1,N4 
280  ?N  (I)  =  PN  (I)  *S 

290  DO  300  I  =  1,  N 

3  00  PW  (I«  (N  3*  1)-:j3)  =  1.0  >  ?'.i  (I  »  (Ni*  1 ) -N3) 

InEVAL  =  -1 

CAi...  BATINV  ,N3,  J1  ,  PUL) 

i.r  (J1.GT.0)  GO  TO  3  50 
oO  TD  440 

310  DO  320  I=1,N 
3  20  3  AV  ^  ,  I )  =  Y  ( 1  ,I ) 

DO  3  4  0  J  =  1 , N 

B  =  CPl^D.IAa  1  (EPS,  DABS  (SAVE  (9  ,  J)  )  ) 


Y  ( 1 , J)  =  Y  ( 1 ,  J)  ♦  ? 
D  =  A  (  1 )  •ri/n 

^AL..  (.,  L  f  J  \  ^ 

DO  jjO  I  =  1,N 

Z  (  N  6 , 1 )  ) 

3  30 

PMI*  (0  -  1)  *N3)  = 

(S  A  V  E  ( ;i  5  1 ,  1 )  -  D  A  V  E  ( N  1  *•  I ,  ! ) 

; 

3  40 

Y  ( I,  G)  =  3A  VE  (  S,  0) 

GO  T U  2^0 

3  50 

12  iBF.NE.O)  GO  TO 
DO  3  0  0  I  =  1  , 

3  70 

360 

3AVC(i,I)  =  Y(2,I 

)  ~  3<i  /  a  ( N  1  ♦  1 ,  1 )  •  :i 
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JO 

10  4 

1 J 

3  70 

DO 

330 

I  = 

1  ,  N 

3  30 

JA  VZ 

{ N  5  4  i 

,  1) 

=  ■{  (2,  I)  -  JA7Z  [li  14I,  1)  *H 

C.-iLL  S 

Ok.  VZ 

(.« 

,?  J  L  ,SAVZ  (N  34  1  ,  1 )  ,2) 

DO  399 

0=1 , 

A 

3  99 

SA  VZ  (9 

,  J)  =.c 

(0) 

4  10 

3T 

=  3 

DD  4  20  I  =  1,'l 

aM)  =  1  (1,I)*-A{1)*3A/2(9,1) 

=  T{2,I)  -  SAVZ19,:) 

Z3tfU2(IJ  =  Z3PCP  (I)  ♦  3AVS(9,I) 

IF  (0A2S  (oA  VS  (9  ,  I)  )  .  LZ.  (iND'-Y.lAI  (I)  )  )  ‘••'1  =  N  I- 1 
4  20  cO.JTINO^ 

I?  (.ir.LZ.O)  GO  TO  4:^0 
4  30  CONTIS  02 
440  T  =  TOLD 

If  t  (H.  LZ.  (  hai'I-»  1.  OOCO  1  )  )  .  A  :<D.  (  (I  WZVAL-MTYP)  .  LI. -1 )  )  uu  T 
IF  (  (ar.  Zw.  0)  .  Jo.  (I«  ZVAL.  3£.  0)  )  3ACL1M  =  r.  AGO  1*  .  2  5D0 
IViZVikL  =  :ii 
I3ET1  =  2 
3U  TO  7‘jO 
4  60  KFLAG  =  -3 
4  70  DO  4d0  X  — 

DO  4  d  0  J  =  1  ,  K 
4  60  =  OAVZ  (J,I) 

K  =  HOLD 
30  =  30OLD 
JSTA3T  =32 
RSTUR3 

4  90  D  =  0-  J 

DO  5  00  I  =  l,.i 

5  00  D  =  J  4  {Z3F03  (I )/YMAX  (i; ) «*2 

I  " tV  AL  =  0 

If  (D.  or.  Z)  GO  TO  540 
IF  (K.LT.J)  JO  TO  520 
DO  5 10  J  =  3,  K 
DO  510  I  =  1  f  3 

510  Y(J/i)  =  I(J,I)  ♦  A  (J)  »E&30d  (I) 

5  20  .-(FLAO  =  4l 
HNE2  =  !I 

IF  (i.DJU5.k.E.  1)  JO  TO  550 
10003  =  LD0J3  -  1 
IF  (  IDOUb.  J  X.  1)  JO  TO  700 
00  5  30  I  =  1 ,  :j 
530  3,k7Et10,I)  =  ZEACE  (I) 

JO  TO  / 0 J 

540  KFLAJ  =  -  2 

I?  (3.  ^Z.  (d.'ilN*  1 . 00001 )  )  Jo  TO  740 
T  =  TOLD 

■1-*  (KFMAJ.h.r.'*5)  jC  To  720 
5  30  ?F2  =  iD/-j  »*Z-;;2*  1,  2 
P.-<3  =  1.^4.,  J 

^  ^  ilfVtvDcfE)  wOci.  ^  'J  ^  TO  d7J 

D  =  0.  0 


0  4  oO 


i 
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DO  5o0  1=  1  ,  N 

560  D  =  J  »  (li:330n(I)  -  SA  7Z  p  0,  t)  ) /'/li  AX  (I)  )  2 
?FJ  =  (D/ZJ  ?)  «*EN  ^3*  1.  li 
570  Pal  =  1.  i  +  ZO 

IF  P'JQ.  L£.  1  }  GO  TO  540 
D  =  0.0 

00  5  80  I  =  1 ,  U 

5  80  D  =  0  ♦  (  X  (K,I ) /7SAX  (I)  )  »  *2 
?21  =  (D/£J«.0  •♦ESC1  *1.  3 

5  90  CONI  IN  J£ 

IF  (?S2.L£.P23)  GO  TO  6  50 
IF  (?R3.Lr.2a1)  GO  TO  660 

6  00  2  =  1.  J/A;1Aa1  (PF  1  ,  1.  E-4) 

N£,»v  =  NO-1 
6 10  IDOuU  =  10 

I?  I  (KFLA^.  £.i.  1)  .AMD.  (F.LT.  p.  1)  j  )  GO  TO  700 
IF  (.iE/i^.LE.N2)  00  TO  630 
DO  o20  I  =  ip.l 

6  20  I  (  NE-•^^«•  1  ,  I)  =  £5  r*OH  (I)  *  A  (K) /DFOOAT  (X) 

630  K  =  NEN^  1 

I?  (KFuAG.  1)  GO  TO  670 
NACJd  =  aAOJO^F 
II'ZI  =  3 
GO  TO  750 

640  1?  {NE«2.  E^.Nv)  GO  TO  2  50 
N  0  -  N  E  W  j 
GO  TO  170 

650  IF  {  ?K2.  Gl.  ?ii  1)  GO  TO  600 
NE2Q  =  Nw 

?.  =  1.  J/A8AX1  (?R2,  1-  E-4) 

GO  TO  olO 

6  60  R  =  1.  J/Ai1  A  X  1  3  ,  1.  E- 4  ) 

NE'rt'iyi  =  N  ♦  1 
GO  TO  610 
670  13  ET  =  2 

F  =  DNIN  1  (S  ,ii;UX/DAS8(H)) 

H  = 

ZNErt  =  H 

IF  (  EO*  N £  «  2)  10  6  80 

HO  =  NEXg 
GO  TO  170 
680  5  1  =  J.  0 

DO  6  JO  J  =  2,  K 
21  =  ai*u 
00  6  ;>  U  I  =  1  ,  M 

6  90  /IJ,1)  = 

looua  =  K 

7  00  DO  7  10  I  =  ],N 

710  I.!AX(I)  =  J.1AX  1  (  t«AX  (I)  ,  DAoJ  (i  (  1  ,  Ih  J 
o'  ZTA  FT  =  N  ^ 
r  ZT’J  it N 

720  17  (NO-Xw-lj  00  '70  7o0 

CAoo  DiFZ'j;.  IT  ,  I,S  AVE  (M2, 1 )  j 
F  =  H/.IOoJ 
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Dj  7  30  :  =  1  ,  ;j 

Y  I  1 , 1)  =  irt  V  O  (1  ,  I) 

JA  /o(2,I)  =  :iGLD*5AV=;  1) 

7  30  :  ;2, 1)  =  3«V0  (2  ,  I)  *3 

•'  0  ~  ^ 

K  F  LA  J  =  1 
00  r  0  17  0 
7  40  KFLAO  =  -1 
KNZrf  =  a 
JSTAuT  = 

FETU  3N 

750  SACO.I  =  OOAXUOABS  (H;1I:I/H0i,J)  ,2  AC'J.1) 
LACO.'I  =  O.IIAMPAC'J.'I,  0A3S(H;-iA;</HGLJ)  ) 
S  1  =  1.0 
Do  7  60  J  =  2  /  K 
n  1  =  H  1  *2  AO U 
DO  7  6  0  I  =  1  ,  'I 

760  7  (J,  I)  =  SAVE  (J  ,1)  1 

li  =  iiOLD-»f.AO'J;i 
70  7  70  I  =  1  ,  N’ 

7  70  1  (  1  ,  11  =  SAVE  (1,1) 

1  ijO  'JO  —  i\ 

00  ro  (130,250,640), ID ET1 
7  30  KFLAO  =  “4 
GO  ro  47o 


SJDEOUIINE  DIFFUN (1, YY ,DTt) 

IdPLICII  .i2AL*3  (A-H,  j-L) 

OXOEl^SUK  0  YY  (23)  ,  YY  (9 , 20 )  ,  A  (20 ,20)  ,  3  (2 0  , 3)  , C  ( 20)  ,  Xi  (  1 )  ,W  A  (3)  , 

1  .iAii(20) 


'  •  ♦  ♦  *  *  *  : 


S012oUTI.i2  CALLED  3Y  CIFSOD  LVAL'JAIZi  THE  DEFIVAIIVES  uF 
DEPELDE.NI  VARIABLES  ST05ED  I.'l  Y  Y  ( 1  , 1 )  ,  I  =  1  ,  2 A  ,  ..3D  STOi-.EE 
THE  OESIVATIVES  IM  THE  ASSAY  DYY.  1  15  IHE  IG  LED  EN  DEG  I 


/  A  X  A  3  u  £i  a 


L  A  T  Z  u  N  A  L  J  X  1 

wN/i-i/coi^ri ,  CO  cf  5 
cj;-;;io;j/l^/co3J2 
V,  3  X  X  0  'X  /  la  i  /  .1  f  X  T 
1  ON  /L  4 /A  ,  o 

C  J .  tk'l  J  N  /  la  7  /  N  J  .k.  J  f  N  «  w  /  4j  ? 
11=2 
aaJ  30 

(Xi  ~  I  i 
30  s.  0  ki  ^  ^  N  J  Zi 


w  A «««« J  Xa  i  Z 


H  —  /  A  aa  U  N  x  .1 .1  S  J  2*  A  va  —  I  ^ ) 


15^ 

(].  1) 

I  r  ;i=  5  Ou 

CT  =  - J.  *Cj'£?  1*  (A  (1  ,  1)  *X.(  (1  )  ‘•A  I  1  /  ♦’••  J  /  A  (1 ,  1)  -  A  i  1,  .,;)  /«  t  1 ,  1;  ’l 

DA  *+D 

C  (I )  =  3(1+1,  1 )  »CT  +  3  ( 1+  1  ,;ir)  *C  Ji.i2 

U5  COiNriNUi 

“ jTASi  I  Hi  DE!<I7ATIVZj  IN  AaShi  D'['i  — - 

DO  1 iO  J  =  1  ,  J 

DTi  (J)  =0  (J  +  1  ,  1  )  *.i  X  ( 1)  +  3(3+  1  ,  ;JT)  *1  - 
DiY  (J  +  X)  =C  iJ) 

DD  110  X=  1  ,  ;i 

D  i  i  (J  )  =  0  1 1  iJ)  +3  (0  +  1  ,  X+  1 )  *ii  {  1 ,  X  ) 

J  iY  (J+j)  =3  YY  (J  +  :>)  -3-  »3  ( J  +  1  ,  1)  =C  0  37  1  *  A  ( 1  ,.<  +  1 )  #  YY  (  1 ,  X  j  /  A  (  1 ,  1  j 

1  +  (3  (J  +  1,  X+l  )-3  lo  +  l  ,  1)  *A  (1  ,S  +  1)/M1  ,  1)  )  »YY  (  1,  .w  +  -0 

110  OJ.ITINJi 

uYY  (J)  =Y0E  ■  1  OCOEES^DY  Y  (Jj 
jYY  (J  +  N)  =  J  Yi  (J+'1)  *COiF5 
120  COi-ii’iN.Ii 
3  IT'J  s;.’ 

EA  D 


j'JOAOOriNE  i'EDEAV  (I,  Y,?X,.'U) 

IMPLICIT  S£AL*3  (A-H,  Q-3) 

OldENSiO:!  1(3,20)  ,?W(400)  ,A(20,2J)  ,3  (20,20) 


^  *  a  : 


««  u  t 


:  »«  a  «  «  a  «  , 


C  EUoAoOTIjL  cAL-ED  3Y  DI?3'J3  LT03E3  THE  PARTIAL  L 
C  L'Z;<I7AiI/EL  jF  THE  DIFFZSENIAL  LlJAlIONj  PPO'/IjSJ  IN  C 
C  oUBROOTIJE  DIFFTM,  THE  PARTIAL  DE3IVATIVES  (THE  JACiaoIA.I)  C 
0  /lEPE  EVALUAIEJ  IN  THE  LA3T  EECTION  IN  CH.vPTER  5- 

J******‘*******»*******»**»*****aa**«aa»**4»<i»**a**=»»»a****»a»«a 


^-OHMON/L  1/COE?  1 ,  CO  EF5 
CuHNON/L J/N  ,  NT 
CC'1.1oN/Lc/A  ,  3 


Do 

1  J 

l\  A  =  1 

,  N 

DO 

4.0 

1m  = 

1  ,N 

PN 

1  II 

+  (2* 

K 

2)  * 

N) 

=C  CL 

F 1  ‘CuiFO «j(ii  +  i  ,a:-;  +  1) 

-> 

b  H 

(II 

*•  i-* 

K  .\- 

1)  - 

N) 

~  c  c 

Ef5  »  J  .  ■‘0  (II  +  1 ,  1  )  »COLr  1»A 

^  1  ,  A.\  +  1  ,  /A  ^  1  ,  ‘  i 

?.» 

(II 

+  i 

^  :i  * 

(-’■ 

■<;< 

-2)  * 

N)=0. 

^  J  pi 

(  II 

+  2  *N 

(2« 

X  K 

-1)- 

N  )  =  ^  D  ( 1 1  +  I  ,  X  X  +  1  ;  -  i  ^  I  I  *  1  , 

1  )  *  A  (  1 ,  X  +  1 ,  /  A  ^ 

10  C‘^  .i  Y  IN  J I 
-'.ET'J  N 
END 


n  Cl 


I 
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subaouri.iz  DiCo:i?  J2) 

DI.l£.iS10.'i  A  v-^0 ,20)  ,UL(2a,20)  ,  3C.\b££  (2  0)  ,  IPO  (20  ) 


:  «  I 


S0330(JTINE  CALi,ED  B'i  aATINV  PEfifCaJS  A  rlSST  STAJE  CF 
GAUSSIAN  ELIrtlNATICN  POUTINE. 


C 


1 

2 

J 

4 

S 


1  0 
1  1 
12 

1  3 
1  4 


1  5 


coa^^oN  IPS 

J2=l 

N=NN 

DO  5  1=1,  N 
IPS  (I]  =  I 

P0  2N  ai1=  0.0 

DO  2  J  =  1  ,  M 

U2  (I,  J)  =  A  (I,  J) 

I?  (SOrt  NPM-A  ES  (UL  (I,  J) )  )  1,2,2 

nOWM&a  =  A3S  (UL(I,J)  ) 

CO  NT  I N  U  E 

IF  (RoiiNaa)  3,4,3 
SCALES  (I)  =  1.0/ROW  NSa 
GO  TO  J 

CALL  SING  (1) 

J2=-  1 

SCALES  (I)  =  0.  0 
CONTINUE 
Nnl  =  N-1 
DU  17  K  =  1  ,  N  a  1 
BIG  =  U.U 
DO  11  I  =  K,>1 
I?  =  IPS  (I) 

SIwE  =  A3S  (aL(I?,K)  )  ^SCALES  (I?) 
IF  (SILE-BIG)  11,  1  1,  10 

B..>j  = 

iJXPIV  =  I 

..ON  ^INUkt 

IF  (DIG)  13,12,13 
CALL  SING  (2) 

J2=-  1 

GO  TO  17 

IF  (IDIPIV-K)  14,15,14 
J  =  I?3{K) 

IPS  (S)  =  IPS  (IDXPIV) 

i.Pi  (I  JIPIV)  =  J 
K  P  =  IPS  (K) 

PIVOT  =  JL  (;<?,;?) 

SP1  =  K*1 
DO  1o  I  =  K?1,M 
IP  =  IPS  (I) 

C..1  =  -JL{rP,  K) /PIVOT 
JL(I.-',K)  =  -F3 
DO  lo  J  =  K?1,N 
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=  yL(i?,j)  (ft?, j) 

lb  COUTINUZ 

17  ^wO  T  iN  J  Ct 

ft  ?  =  I  ?  i)  I  .'1 ) 

17  (JL(.<?,Mi)  19,18,19 

18  CALL 
J2=-  1 

19  PETUBN 
END 


GUuHOUTINE  SOLVE  (  NN  ,  PUL,  3, a) 
iNPLICil  ft LAL*8 (A-H,v-2) 

DI.l£NSx3;j  ?UL(20,20),  3  (2J),a(2Jj,  IPS  (27) 


C  SU3&  Oil  TINE  CALLED  3  2  OlrSUB  ?E1<F0B:1S  THE  SECOND  STAOE  C 

C  (THE  BACft  SUESrirUTICtl  PROCESS)  0?  THE  JAUSSIAN  C 

C  ELIHINAIIUiN  EOvJTINZ. 


t  v  **  **  m-f  I 


COaHON  iPS 
N  =  NN 
N?1  =  N*1 
IP  =  IPS  ^  1) 

X(l)  a  a(I?) 

DO  2  I  =  2 ,  N 
I  P  =  IPS  (i) 

I H  1  =  i  -  1 
S  U  H  =  0.0 
DO  1  J  =  1,I.'11 

1  SUE  =  Sd.l  ♦  PUL  (IP,J)»X(J) 

2  ;i(i)  =  8(iP)-su;i 
I?  =  IPS  (N) 

.<(N)  =  X(..)  /PUL  (IP, N) 

DO  U  I  SACK  =  2,N 
I  =  NP I'idACft 
I?  =  IPS  (I) 

I?1  =  !♦ 1 
S Ja  =  0, 0 
DO  J  J  =  I?1,N 

3  S  JH  =  SU.l  +  PUL  ( IP,  J)  *.<  (J) 

4  X(I)  =  (X(x) -SUN) /?rjL(I?,I) 

psru  dn 


SUoaoUriNE  sing  (IftSY) 


'  ««  «4i«  «  < 


:««««««*«! 


I  •  *  a  a  ' 


SOSrOJTI.it  CA-LZD  JY  DZCCNP  iNOICAiES  THE  ZaSOaS  IN 


i 


Ll  O 
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rji  SA'JJSIAii  ELIZ-INATIOW  aOJTINZ. 


c 


1'  FC/ jj  Ai  ( 5  4  .i  J  .'i  .\Tr>  I X  WITH  ijOii  ill  3ZCUMPO.TE. 

12  r  0  A;!  AT  (  D4  ri  3  SZ  N  J  HL  A  F.  ilATSIX  U  DECCiiPOii  2.  ZErtO  DIV’IZZ  13  _iCL7E. 

:  i;j  IKPAUV.  ;iArFix  u  nzafly  z^ajela:-. 


13 

FOLHAr 

l54  H  J 

30  CONVERSE 

NOUX  = 

3 

SO  ro 

(K2, 

3)  , 

IWH  Y 

1 

wair  2 

( N  J  u  r 

,11) 

so  TO 

10 

2 

WRITE 

(300I 

.12) 

o 0  TO 

10 

3 

wall  E 

t  u  u  r 

.13) 

10 

LExuau 

ill  ^  0 

s'jEaou 

:iuE 

HAT  I 

:;v  (TN.N 

Di.lEUS 

I  J.1  ? 

W (400)  , A  (20 

««  «  «  * 

3  U  3  j  J  J  T I 

.IE  w .A L L Z D 

3Y  DIF 

1  •  *  »c 


L 

c 


GAJoSIAN  2  L 1  i'l  1  i'l  AT  I  0  N  liOUTIliE  OF  T;l£  ^lATFlX  jT'JAZO  1 3 
TH2  AS  FAT  Pa. 


DO  1  I  =  1 ,  :i 
DO  1  J  =  1 ,  ;i 

1  \  (I,  J)  =  Pa  ii4  {j-  1 )  *:;3) 
CALL  DEC^/*P  (*4/AfUL«J1) 
n  ZTU  RN 
END 


jUwoJUlI.lL  CUFZNT  (CD,A,Y) 
1:1.^1011  R  jAL*3  (A-K,  0- Z) 
DI.'JE.'iSIj.V  A  (2U,20)  ,Y  (8,20)  ,  Xi  (2) 


C  SrJDSJUTIJZ  CALCULATES  THE  C'JSaEUT  JEUSIl'i  31  EcUArioL 
C  (5o)  . 

c  c  o=caa  hzu:  jl  ..o  i : y  ,  Ar  ?/c  **:  : 

C  A=  DIoC  L  i  I  A  i  lu  N  C(3S?FICIZNr  .''.ATZiI.C,  A.  C 

u  AT  C0I*««0CAr.A.\^.4  201. 1 73*  ^ 


n  r  .  n  n  r,  n  n 
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CO  j  j/  T 

Co. -I.-l  ON/ /  \0 
Cj;iAON/„o/C  Ju?-. 

co=-cjc;:-’4»  ( ti ,  1)  *X5  (1)  ♦■A  { 1 Z)  ( 
Du  3  0  J  =  1 ,  ;i 

CO=  Cu-CUtZ'  4«  A  (  1 ,  J  4-1  )  Z  (1  /  J) 

50  COMTINOZ 
z  ciru  till 
END 


o  O'  3  S  O  Li  r  I  .i  E  C  0  ‘.J  C  r.  N  ( C  1  ,  C  2  »  N  \  ,  U  D  ,  D  ZT  0  0  1 , 0 1  ?  1  /  : ) 

I.'IFLICI:  SZAL«3  (A-H,  0-C) 

DI.I'ZNZ.D.;  Cl  I  60)  ,C2(60)  ,P.OJT  (/O)  ,  Jl  J  1  (2  0)  ,  'Z  (  1,  .:.0)  ,  E3  NZ  Z>  (2U) 


U 

C 


SOSf.OJ  I’x.j  Z  ZVALOAIZS  TOE  CO  .iC  Z.or  a.U'I  0  N  AT  ZAC.'i  DZCIZZC 
IN  X  OINZCTION  E'£  AN  IN  I ZN  ?0  OA  Z 10  N  SOOTI.iZ  OCINC 
THE  KNOO.i  ZJv<CZNTRATIONS  AT  THE  COOLOCAIIOl  POZNZS. 

“““  i'Ar.iV.lii.XiR.j  — 

C  1  ,  C2=  CD  NZ  E.i  1  Zii\  T I  ON  S  u?  N'lTRAIZ  AND  .n'ZDRJX'Z  lUNJ  , 

N  Eo  2  ZCI  1  /  ZL 1  /  IN  D  I  HE  NO  lONL  Zoo  IN' IT. 

DETA  =  a  INXZqZA-,  OlilENoIOJLZoS,  0  <  DZT.-i  <  1. 


L 

c 


c 


CO:i.10N/L2/COZZ2 
CjCjiaN/LO/N  ,  N  : 

C0'13uN/L3/Xo 
DO  3  0  J=1  ,;i.( 

ZIA  =0iTA’‘  (J-  1) 

u  A  L  I N  i'  o  D  ( .<  D  ^  N"  T  ,  ET  .A ,  R  GO  T  /  D  X  c  1  ^  .0  N I  ^ ) 
C)  ( J)  =iiN  Z2  (  1)  *X5  (1  )  4XINT?  (NT)  -1. 

C2(  J)  =.XiN:2  (  ))  *Xo  (2)  ^XIN  :?  (NT)  ‘COZZO 
DO  7  0  K=  )  ,  .< 

-  1  (J)  =C  1  (C)  +XI'M?  (K*  1)  *l  (1  ,K) 

7  0  Cx  (J)  =02  (o)  ♦•XINT?  (X+  t)  (1,  t. 

30  Cu  .i  r  IN  u  z 
NZIUaN 
ZND 


T  or  7«'J A ,  0  o  i  .'i  ai  I ; 

I.l  2'-  IC  12  J  Z  .i«  ♦  o  ( A  -  li,  0-  Z ) 


O'J.  .^UwaX.tM  X.(^.^^■'i^.aO  %d  •J  f  .(.<0  .  .lu  Oa.iaaO 

a'  A  .k  .'tx  a  o  aO  ^  a.  .3  A  L  L  li.A  ak  ji:*  lO  J  o  a  0  Z  H  ii#  x  .<  Co  x  A  o  ..  a  .i  a^ 

oPE^Ow  201.'  i* 
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^.\'!Z6=::.6 

IF  (u6-  jZaJ*  Oj  1 .  .\ND.ZS>  LT<  J.  'jlJ) 

If  (-J..JG.  J-  JI  O.AND.ZS-  LT-  0.  100) 

Ir'(u3*vjZ«0«  1j0)  Du  S~  0*0  5 

uS  —  uS^OOu 

FACro»l=^0/SAVD3 

CO  ZF  J  =  OD  if  j/r  ACTOR 

CjZF-4  =  CJiF4/r  ACT0  3 

CO  if  D  =  OJZr  5/F  ACTOa/?  ACTO.' 

DTAJ  =  DrA  J/ 1 0 j . 

Z:\J.iL=.  00  1  *JTA7 

DlAJ.1A  =  2.0*Di  A  J 

JiTAaT^O 

FE-U3N 

Fl.D 


OijO'^O  •  0  005 
DOS=0 .005 


SUoaOCriNZ  ZXPAA'D  (I,  foot,  Olr'l  ,f  ACTOa,  NO) 

I.lFLICir  .•vZAu*5  (A-H,  0-0) 

DI.-1Z..5.0-,  i  io,  20)  ,  ROOT  ( 20  )  ,  DI  Tl  (20 )  ,  ^  i  T  ?  ,  2  0 ;  ,  0  I  ( ^  C  )  ,  iO  1 2 ) 


'««««*«««« 


t«*#«»***^i**»i***;" 


S'JH30UTI.4Z  EVA;.0ATES  THE  CONCSNrSATIONO  AT  TNE  COLuCCAIlUA 
P0I:J7S  in  iiiz  COOKDIKATZ  DOE  TO  THE  CHANGE  Of  THE 

oPLiNE  POINT.  C 

C 


COH;iON/L2/OOiF2 
C 0 HO O N / L J /  N  ,  NT 
COHNON/La/iCO 
00  7  J  ;  =  1  ,  N 
HI  (I)  =1(1,1) 

7  0  i  I  ( 1  ♦  N  )  =  ■£  (  1  ,  I  N ) 

DU  100  1=  1  ,  N 
RE=aOJr  (!♦  1)  *FACTC3 
IF ( HE . GT. 1 . 0 )  00  TO  50 
C.il.1  INlfi'  (ND,  NT,  3E,  .1001,  OIF  1  ,  TIN  TP) 

I  (1  , 1/  =  1 1:.  .r-  {  1  )  *XS(1)  <■  il  ..TP  (..T)  ■'I.O 
I  ( 1  ,I»N)  =  ;i.iT7(  1 )  ♦X 3  (2)  >  YI  ..TP  (Nl)  *CJEF2 
Du  ”  J  r.  =  1  , 

;'  (  1  ,  .  i  =  i  I  1  ,  D  ♦’TINT?  (!<♦’)  £  ( H) 

60  I  ^  1  , 1  ♦■  (i )  =  4  {  1 , 1  ►  N)  1 1  N  r  ?  V  H  *■  1 )  *  0 1  ( IX  N  ) 

Go  .0  100 
9  0  :  ( 1  , 1 )  =  1  .  0 

4  I  1  , 1  N )  =  C  0  i  F  2 
100  COnTINOE 
lElOKN 

M  .<  U 


n  n 
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Fj;iOTioN  AjAi(iC;(,K,iArt) 

I.'iFLIOir  5  (A-H,  i-O) 

DIlEiiSiOM  <\  (^0,20)  ,3  (2  0  ,2  dj  ,Xa(1)  ,2Ari  (20) 


c 

c 


««««««« 


C  AOXI  IS  Th£  1<AJS  OF  THE  FUNCTION  CALLED  31  ZSYST.1  IN  C 
C  THE  SUdfiJUIINa  DIFFUN  TO  F'JSwISa  THE  VALUES  OF  Ej'JATUNS  C 
C  (52)  AND  (53)  .  C 


I  «  < 


CC.13ON/Ll/2  0a?  1  ,C0EF5 
COMHON/LZ/COEF2 
C0;i.1UN/LJ/N  ,  NT 
COili'lON/L'i/A,  3 

COi'lHJN/Ld/— o£F3,C0  N1  ,C0N2  fSA.lAfBSTA 

;A3  (2U)  =-3.  OCCEFI  *  (A  (1  ,  1)  *XX  (1)  ♦  A  (  1 ,  NT)  *  1. 0)  /A  (  1  ,  1  y  -^  (  1  ,  ;.T)  /  A  (  )  ,  1  j 
1  *CJEF2 

DO  5  00  J=  1  ,  N 

lAa  (20)  =  w  AN  (20)  -  3.  ♦COE?  1*  A  ( 1  ,  J+  1 )  •2A3  (J;  /  A  ( 1  ,  1)  -  A  ( 1 ,  J<-  1) 

1  {J  +  N)/A(1,1) 

5  00  CONTINUE 

IF  (^ Ao  (20)  .  LT.  3.  )  CA?(20)=0. 

AUX1  =A  (  1  ,  1)  *  Xi  ( 1)  ♦  A  (  1y  NT)  *1  .  *•  (CO:!2  *  (vAE  (20)  »«  UABA)  -CON  1  * 

1  (iX  (1)  •♦SETA)  )/COSF3 

DO  t)00 

600  AUi1  =  AUX1  trt  (  1,  J+ 1)  •QA5(J) 

FETJ3N 

END 


FUNCTION  A  JX2  (:<:(,K  ,3AS) 

IkiPuIwir  EaA^*3  (A'~rif  0'*'^) 

DIMENSION  A  (<:0,20)  ,5  (2  0,2  0)  ,  I  (3,2  0)  ,  X  X  ( 1 )  ,I;Aa  (  1) 


'ii«  *#«  «  *•« 


C 

C  A 0X2  IS  THE  NA.iE  OF  THE  FUNCTION  CALLED  oT  ESYETi  Ii« 

C  THE  EAIN  PEuGSAH  TO  FURNISH  THE  VALUES  OF  Ei^UATIjNS 

C  (52)  AND  (53)  - 


C 

C 


CO. -I/.  ON /i,  1/COEF  1  ,CO  EF5 
C  0  EH  0  N/  L2  /C  U  E  F2 
C  jii10N/L  J/N  ,  .JT 
0  0  O.'i  O  N  /  L  4  /  A  ,  d 

Cu.1.1  ON  /  L  o/  -  0  EF  J  ,  CO  N  1  ,CO  N2 , 0  ACA  ,  oET  .A 
Cj.,.10N'/LU/i 

3na(  1)  =- J.  *C.,EF  1*  (A  (  1  ,  1  )  ‘Xa  il  j  ♦  A  (  1  ,.4:)  *1.0)  /  A  (  1  ,  1)  -  A  (  1  ,  NT)  /A  l1  ,  '  ) 
1  •CuZF2 

.30  5  00  J  =  1,N 
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3A5  (  1)  -5A.^  ;  1)  -  3.  *COE?  I*.;  i  1 ,  J  +  1  j  *  ,  J)  /.I  ( 1 ,  1  j  -  A  ^  1  1 , 

5  00  cooin.iiis 

IF  ll)  .  LZ.  0.)  3.'F[l)-0. 

A  U  a2  =  A  (  1 ,  1 )  *  A  A  ( 1 )  4-  A  (  1 ,  T)  « 1  -  (CC.i2  *  (  3  AA  ( 1 )  *  *  Jn  .lA)  -CJ;.  1  * 

1  U-U  '■)  **3FTA)  }/COEF3 

DO  6  JO  J=1  ,  .< 

6  00  AUi2  =  A  Ji2  ♦A  (  1 ,  J  +  1)  *Y  { 1 ,  J) 

5  JT’J  PH 
EHO 


10<i 


'  C 

c  rnrs  ?iijo:tA/  juL.7^5  j a t  10 it j  ( 4 d ) ,  ( ^ 7 )  ,  and  {^2)  1  r. d 0 j •-» n  w 

C  (55)  «i.Ta  i«j  INITIAL  CCNDITiONS  JJINii  SPLINE  COLLCCAriOS  C 

c  r.srHOD  t'ja  i;ie  case  cf  the  zlectsolite  containins  cadhiu.i  c 

C  IONS.  IN  THIS  CASE,  T;1S  TWO  SSaAOUriNES,  DIFF'JN  AND  PEDaSV,  C 

C  CALLED  BY  THE  INTEGRATION  saaBOGTlNS,  DIFSUa,  HAVE  TC  BE  C 

C  CHANGED.  HEAN WHILE,  THE  3AIN  PEOGEAH  AS  WELL  AS  THE  CTnatt  C 

C  S GBROaTINES  STAY  THE  SANE,  THE  CONCEN 12  A TION  PROFILE  UF  C 

C  THE  CADHIUH  ION  HAY  BE  OBTAINED  BY  THE  ELECTRONENTPALITY  C 

C  CONDITIO.,.  THE  COBS  ENT  DENSIIY  DATA  A2E  CALCULATED  IN  IHe,  C 

C  SAHE  WAY  BY  S^idATION  (56).  THE  H  E  TLa  JGE  N  EOUS  SEACTION  BATE  C 

C  CONSTANTS  OBTAINED  BY  LAST  P2CG2AH  ABE  USED  TO  IDENTIFY 


C  THE  riO.IOSENEOUS  REACTION  RATE  CONSTANT  BY  FITTING  THE  C 
C  CALCULATED  CURRENT  DENSITY  DATA  WIlH  THE  £a? ERIHEN TAL  C 
C  DATA.  C 
C  C 
-  -  ADDITIONAL  NOHEHCLATURE  -  C 


C  SH  ALLK  =riO;‘lO  SB  NEOUS  REACTION  RATE  CONSTANT,  1/3EC.  C 

C  DaiaK=0lHENSlON  LESS  HO .10 G EN EOUS  REACTION  RATE  CONSTANT,  C 

C  0£  i  0  i\  =  S  Hrt  LLK  ♦  D  .‘In  X *  •  2/5.  2b  E ~5  .  I 

C  DKS?=Di;iaNSlONLSSS  SOLUBILITY  PRODUCT  0?  CADHIJ,'.  HYDROXIDE,  C 
C  DKSP  =  2. Oa-23/C1 a**i.  C 


V. 


laPLICIT  REAL*8  (A-H,  C-2) 

DIMENSION  A  (2  0,2  0)  ,S  (20,20)  ,C)  (60)  ,C2  { faO)  ,  C  (23)  ,  Y  (o,20)  , 

1  0  IF  1  (20)  ,DIF2  (20)  ,DrF3  (2  0)  ,  ROOT  (20)  ,  ‘/ECT  (20)  ,  ZnR-R 

2  Y  HAX  (20)  ,S  AV  E  (24, 20)  ,P.i'  (400)  ,  X3  (2)  ,  XX  (  1)  ,  4n  (  J)  ,S  AR 
EXTE  RN  AL  a  J  X2 

COS;iON/L1  /wOEFI  ,  CO  EF5 
CC  '..ION  /L.:/COi.F2 
C  *!10i.  ,j  ,N. 

C  .aX  vj.)  /  L  4 /  n  ,  L 

COilHON/LS/XS 


C0:lH0N/Lo/C0E?4 
CO  .'l.'l  ON  /L  7  /  N  Si  G  ,  N  Z2  ,  ZP 
CU.H.’10N/Ld/COEP3  ,CO  NT  ,CO  N2  ,  S  AHA  ,  BE  T  A 
COKHON/LR/. 

CUH.ION/L  10/0dIGK  ,  0  KS  ? 


INPUT  DATA  ” 


n  an  j  ( 5  ,  J  'i  0 ; 
-R  SAD  (5,320) 
R  a,'»  D  I  5 , 3  3  3 ) 
P  BAD  i5 , 3  4  u ) 


D,  j,nO,N1  ,Ai.,aM,uPS, 

,.1  3,  C2B,  D;1  AX,  DTISS,  NY,  SS 
.j  A  .1 A  ,  3  c.  T  jA  ,  C  0  N  2  ,  C  0  n  1  ,  S  .1 L  L  H 
U  3  ±\j  f  N  c.^  ,  E  P 


•  .ALC’JL.ATE  SaHa  CO.iSTANT. 


:  J 


?  R  a  S  5  n  .1  - 


CDEF1=  (1.  :>02D-5)  /5.  26D-5 
ZOZe  2-CZs/C  1  3 


n  r,  n  ri 


r 
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Cwif  1=  t  1 . 's  J^3-5)  *  Cl  2/D.‘.:\  X/ZJ 
LJZFU=2.*  I  U  'J  )2D-5)  »C13*'>T4cj7-/J.'i.‘iX/2J 
—  *5— jS 

3 3 4  j  j  44 rt  C  L  X *  l) 4X4^ .C  *  A.’.  /  5*  2 o  D 
DKj?=(Z.0-2  jj  /{Cia*»3.  ) 

2X=241Aa/(4'Ia-1) 

3:AJ=  (5.  4Li3i-5)  *DTI 

TAC'1  A.X=  (5.  2o0-5)  »T  i:i  A  2/D  3  A  .( 

5JTA  U=DTAU 
02rA=i./v:i.c-i) 

M=N  ♦■;(  0  +31 
:i  2A=  2*4i 

-  JZr  INITIAL  0I.12N  TIAZ-'J,  - 

TAU= J. 


J  4  4.'  .i.  4t  4J  i. :  1  w 


^lAxi^ua  AND  4ii;i::iJ:i  dlxendlonlldd 
Il4'1>^  I  .<  vZ  3  H  4IE  NT  To  ’.JSED  til  THE  IN  TED  ?  .\  I  lo  4« 
EUuSoaTINE  OIFSOa  - 


D:aJ,U=  -  0  0  1  'DTA  J 
DT  AO  .1A  -  2  .  J  '*  .\  'J 

•*  A  IT  E  (  o  ,  i  b  0  ) 


-  EVrtLOArZ  THE  COLLOCATIOj  ?Oi;i  IS  (.-vCOTO)  07  Z  AZ 

JAC05I  POLTNOKIAL  0?  OBDE3  0,  AS  WELL  AS  TEE 
Lib  ST  AND  SECOND  DEhiVAIIVES  Of  ThE  P'jLi'NoEIA- 
A I  IriS  POOTS.  - 

•  OEsJb^  (  4^0^  4^1  .1O  ft*  1f  f  OwfDXr  1/  D  x  1 2  f  OX*  4.) 

/•4SirE  (b,J50)  (?. OOi{J)  ,0=1  ,.tT) 


ALCJLATE  TiiE  DISCS  £TI2Ari04'l  CZZEf  FICIE4iT  MAT! 


ID 

20 


10=1 

00  4.0  1=1,4.  I 

CA*.L  J  :  0  2  Z  ^N0,N,N3,N'1  ,I  ,iD,0I?l  ,DIT’2,Dif-l,SCOT,V..,cr) 

DO  10  0  =  1  , N  T 
A(I,J)  =V.CT(J) 

w  u  N  T 1 4^0  Zi 

**“‘"**“‘“  V.4'i4.0oL4'l«W  4..C*  Dx^4«4^WbX44.'(.^^4l  CC4Z44trl4.L.J4lb  .'4.\4  4'.  4.  4  4  0  - 


3  0 

-.0 


DC  *=!,.<  4. 

C  ft  4.  4*  0  4'  .t  4.  4,  ,  4.  ,  t*tJ  ,4.1  f  ^  f  40,44X4'  I  ,Ol^  2,  4/44'  0,4. 

D../  JO  ij  —  (,.1. 

:  (1,0/  =7ECT  (J) 

.J4,  t  T  I.t  j  .i 


I t) 


1 4*  4^  J .  4 E  VA  LO  —  S  jc  .\S'.j>4l4:,.T.>  ^  jE  JSE. 

X  It  4  -  4j  P  \T  I  '..i  N  3  J  j  ZO  J  Z  I .'*  Z  Ox  4*  S  0  3 


i 
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.'1r  =  l 

sJSTA  I'T  —  O 
I5TJ  ?  =  0 
:'.AXDE2  =  7 
DO  47 

47  y.'lAX  il)  =1- -7 
1  =  0 

C  -  iiiDOI  TH2  IHLTIXL  JCNDiriONS - 

C 

DO  50  j=1,:i 
7  11  ,v])  =  1. 

50  Y  (1  ,  J^’J)  =0  007  2 
XS  (1  )  =  1. 

XS  (2) =OOSF2 
WrilTE  (b  ,  J^O  ) 

60  TIO£=TAO*D.'.Aa*DOAX/ (5.  2bD-5) 

V-  ———————  — ALCULATZ  AMD  PRXMT  X^iS  Cu'Ru— MT  DENSITY  AT  ihXj  ri-.'i  — 

J  PSINT  THE  CONCE.VTaATXOiVS  0?  NITRATE  AND  HYDaOXY  10  US 

C  AT  TOE  NT  COLLOCATXCM  POINTS  - 


Call  c o e  s n  x  ( c d , a , y ) 

«i3ITE(6,jo0)  TIME,CD,X3(1)  ,XS(2)  ,  KFi.  AG,  ITEEM  ,  AS 
NEIT£(6,700j  ( Y  { 1  ,L)  ,  L=1  ,  NEE) 

C  -  EVALUATE  AND  PRINT  THE  CONCEN T E A TI C NS  0? 

C  clIi'EATE  AND  HY030XY  IONS  AT  EACH  POINT 

C  IN  X  DIRECTION  AT  THIS  IlrtS - 

C 

CALL  CONCEN  1,0, C2,  NX,ND,OErA,  ROOT,  DIF  )  ,  Y) 

;»I!ITS  (6 ,2  JO)  (C  1  ( J)  ,J=  1,  NX) 

«nITE(o,250)  (C  2  ( 0)  ,  0=  1 ,  NX) 

C 

C  - CHECK  IF  THE  INTEGEATION  TIHE  AT  THIS  .lU.IENX  REACH -3 

C  THE  HAXaiUa  TINE  TO  3E  INTEGRATED  - 

C 

IF  IT  AU.  GE,  TAU  MAX)  GO  TO  190 
35  1  =  1+  1 

I?  (I-  EQ.  1 )  GO  TO  9  0 
IF  (ISTDP. E^.  1)  GO  TO  90 
C 


r* 

c 

c 


TESr  IF  THE  CONCENT  RATI  ON  0?  NITRATE  ION  AT  THE  NTH 
EEuLGCATION  POINT, IE,  THE  LAST  ONE  BEFORE  THE  SPLINE 
PC*.ai'  13  A  IT  HI  N  J.OOCl  01.''.  EN3  I  JN  LESS  C  J  NCEN  IR  .t  » I  u  ,< 
UNil  OF  THE  boundary  CONDITION.  - 


C 


c 

C 


IF  1 1  1* -Y  ( 1  ,  3)  )  .  LT.  0.  OCOl)  GO  TO  9  0 


if  THE  TEST  FAILS,  THE  3P..i:iE  POINT  IS  I  NCI- iAS  Ej  ,  I E, 
IDVi*^  7UHTKE3  INTO  THE  SOLUTION.  - 


V. 

c 
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^  3  f  D  /  r  Av«TOB  /  DTA'J  /  D  ^AJ.IX  ^  D  lA  JXji/ 

J3TA3T) 


■~~~  cVAXidATE  THE  CO  NC  ij.«  13  AT  ION  J  J?  .JITPATE  A.'lj  d/DiiOX'/ 

IJHb  AT  THE  CO  L  LOC  <&  rX  0  N  POxNT3  AT  THE  .'ic/i  CCudJLJal^ 
OJE  TO  THE  i:iC3ZA3E  OF  THE  SPLINE  POINT. - 


CALL  EXPANl)(Y,300T,  DI  F  1,  FACTOR  ,  ND) 


- INTEGPATS  EQ'JATIOMS  (NH)  AN3  (43)  USING  THE  C  J  NC  E  N  13  AT  1 0  NJ 

AT  InZ  COLLOCATION  POINTS  AT  THIS  TINE  TO  OLTAxN  T  .iE 
CONCENTP AT  IONS  AT  THE  SUCCESSIVE  TIHE  - 


90  CALL  DIF  SUB  (NER  ,TA  U,  Y,  S  AVE  ,  OTA  J  ,  D  T  AUH  I  ,  OTA  UNA  ,  EPS,  NT,  YNA.C  , 

1  c>RRGR,KFLAG,JSXARr,NAX0E3,Pii) 

IF(  KFLAS.  LT- 0)  GO  TO  190 
DTA  Ufll  =  .  00  1  •DTAO 
OTAUNA  =  2.  a*DTAU 
I?(I.G£.2)  Gu  TO  93 

IF  (Y  (1  ,  1  )  .  S  1,  ,  OR.  Y  (1  ,  2)  .  ;T.  1..  or.  Y  (  1 ,  3)  .  GT.  1.  .0.1.  1(1  ,  •*)  -  OT.  1  . 

1  .OS.  DAoS  (  Y  (  1 , 5)  -  1.  )  .  GT.  0.  0004)  GO  TO  90 
30  TO  95 
93  N1=:/IC 

a2=ai*ic 

I?(I.Sw..'l2.0S.TAU.GE.TAU.NAX)  GO  TO  95 
GO  TO  d5 

- CALCULATE  THE  SURFACE  CONCENTS ATIOl.S - 

95  OLDXS2  =  XS(2) 

XX(1)=Y(1,  1) 

ITbS  N=500 

CALL  ESYSIU  (aUX2, EP, NSIG, NEZ, XX ,irEHN, 4A,oAR, lER) 

XS(1)  =XX  (1) 

XS(2)  =-J.  *COE?1*  (A  ( 1,  1)  *XX(  1)  ^-A  (  1,  NT)  *  1 . 0 )  /A  v  1 ,  1 )  -  A  (  1  ,  N  T)  /A  ( 1  , 1 ) 
1  ’COEFE 

00  150  J  =  1  ,  N 

XS  (2)  =  XS  (^)  -  3,  *CCEF1  *A(  1 ,  Jt  1)  *1  (  1  ,  J)  /  A  ( 1  ,  1)  -  A  (  1 ,0^  1) 

1  *i  (1,  J»N)/A(1,  1) 

150  CONTINUE 

If  (XS  (2)  .Lr.OL3X32)  IS! 0?=  1 
GO  TO  oO 

190  NRir  S  (a  ,2'40  )  0  (,  3  OIA  U  ,  N  X ,  0  TAJ  ,  D  21 H E ,  ON  AX ,  AL ,  B  E  ,  EG  ,  E  PS  F Ea  j  , 

1  E?,NSIG 

'iiN4Lii{a,245)  b  ET  A  ,  GA  N  A  ,  CO  N  1  ,  C  J  N  2  ,  S  N  A  L  L  X 

- T03.UTS  TON  INPUT  ANO  UJTPJT  STATEXEN  TS - 

230  F0oaAr(/,23X,'N0J-',lX,:»D11.4/3Jx,i011.4/30X,sJi1.'4/j0X,90n.^ 

1  /jO  A  ,  90  1  1  .  4/3  Ox  ,  9D  1  1  .4/3  OX,  901  1 . 4/ J  0  i ,  9  ^  1  1.  4/ jO  .v  ,  J D  1 1  .  4 


J 


A 
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2 

2  i.0 

1 

2 

3 

a 

5 

2  45 

1 

2 

3 

2  30 

1 

2 

3  10*" 
320 
330 

3  40 

2  50 
360 

3  70 
300 
3  00 
7  00 


/JO.i  ,  OD  n  .  4/3 0. i,  OD  1  1 . 4/3  02,  30  11 .4/3  02, 9 D  n,  4/ 30  X,  Ou  1  U  4) 
FOOMAi  (//,5X, '2  I:;TEPVAL=' ,  012.  4,322,  '  INITIAL  lAO  iI52> 

,0  12.4/52,  '  MO  Or'  POINTS  IN  2  OIS=  ',13,202, 

•LAST  TAT  STEP  SI..E=  '  ,  0  1  2.  4/52  ,  '  3  E  AL  11:12  13122732=  ' 

,  0  10.4,242  ,  'DIFOSIOli  LEliOTH=  '  ,  0  1  0.  4//5  2  ,  '  AL=  ',.-7.2,3  02 
•4E=  ' ,F7.  2,252  ,' 2J=  ' , F7 . 4//5 2 , ' E?S=  ',010.2,27a, 

'KFLAO=  '  ,  I3//52,  •  E?  =  '  ,  J  1  0 . 2 , 27  X , '  :i  3 10=  ',13) 

FO:-i;i  AT  (//,  5  2,  'EEACTIOM  OSiOSE  OF  NO]-  ION  =',F7. 2,25a, 

'dSACTION  050E5  OF  OH-  UN  =  '  ,  F  7 . 2//5  2 ,  '  13T  C:)NSTA3T  =' 
,012,4,252,  '2ND  CONSTANT  =',012.4 
//52, ' SEACTICN  RATE  CONSTANr=  ',012.4) 

FOSa AT l/,252,  'OH- '  ,2X,  9  01  1.  4/3  0  2,  901  1 . 4/30X, 90  1  1 .4/302  , JO  1  1 . 4 
/JO A, JO  1 1. 4/3 Oa, so  1  1.4/3  0  2, SOI  1. 4/3  0  2, 9  0  1  1.  4/ jO 2 , 9 D  11 . 4 
/3  02,  JDl  1.  4/3  OX  ,  90  1  1. 4/302,  901  1. 4/302,  901  1. 4/30  2,  ?L  11-4) 
FOS;iAT  (41 3  ,  2/5-  1,  2 CIO.  2,1  5,  .-a.  4) 


F0R;1  AT  (40  1  1  .  4,I7,F7.4) 

FURHAI  k^F7.  2,3012.  4) 

FQR.1  AT  (215,  012-4) 

FOSHAT  (152,  1  OF  10- 6) 

FOR.I  AT  (//,  u  U.  4, 1  2  ,D  12-  4,  22 ,20  3  0.  10, 215,  F2  0.  7) 
FOSH  AT  (/,  2  2,  oF2Q.  6/2X,  6F2a.  o) 

Fo:<a  AT  (22,  '  COLLCCAIION  POINTS:  ') 

.='OR.1  AT  (//,  J  A,  '  Ti;iE  •  ,  12;C  ,'C0.-iSENT'  ,452,  '  lOM  CON 
FORHAT  (302,  10  F10.  5/3 OX,  10F10.  5) 

STOP 


ENTRAriON') 


END 


SUSaoUriNE  OIFFON  (I,  YY,0’2X) 

IHPLISLI  RZAl«>3  (A-H,c-C) 

OlaENSION  0  22  (20)  ,  YY  (d,  20)  ,  A(20,20)  ,  3  ( 20 ,  .;0 )  ,  C  ( 20)  ,  22  (  1 )  ,  «A  (3)  , 
1  ^AR  (20) 


C  Su3R00TI.U  CALLED  BY  DIFSaS  EVALUAIES  THE  DERIVATIVES  OF 
C  OEPEnOZNT  VARIABLES  ST03E0  IN  Y Y  ( 1  ,  1)  , I  =  1 , 2* N ,  AND  SToFES 
C  THE  DELIVAriVES  IN  THE  ARRAY  Dii.  T  IS  TilE  iNOEPE;iOENT 
C  VARIABLE. 


E2TERNdL  AJ21 
COH.ION/Ll/COEF  1  ,COEF5 
CuaHON/L2/CO£F2 
C0.1.10N/L  J/d  ,NT 
C  J  I’iH  ON  /  L 4  /  A  ,  3 
CO Hi'lUiJ/ij7/.*  S.^,  Nr,C  ,w? 
CO/.;u;(/L10/Di5lGK,  OKS? 
11  =  2  *11 

OU  3  0  1=  1  , 1 1 
.-AR  (I)  =YY  (1,1) 

30  COl'ITINJE 


c 
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calculate  T:1Z  7AI.JE5 


j  LI  .1 4  I  i  C  w 


l-<  =  0) 


u.l 


(1/  Ij 

i'rzi<;'j=3oO 

^AvJaI  f  as  A  A  f 

CT=-3.  *«-JEf  1*  (A  (1  ,  1)  *XX  11)  «-A  (1  ,.Jr)  *1.  )/A  0 . 1)  -  3  1 1#:-!)  O  ,  1)  *CCEFZ 
DO  4  3  i=1,;i 

C(I)=u(I*l,l)  ‘CT^Bd^l  ,Mi)  *C0cr3 
45  CO^'jriNJE 

- SXJjE  I'HE  DEEIVATIVE3  Ij  AAHAX  DYY - 

DC  1  2G  J=  1  ,  'A 

DYY  (J)  =o  ( J  *■  1  ,  1  )  *  X  X  ( 1)  ♦d(  J*  1  ,  .iTj  1  . 

DYY  (J*:j)  =C  (J) 

DJ  110  K=  1 , ;; 

DYY  (J)  =DYY  (J)  +3  .<f  1)  *YY  (  l.K) 

110  Q  Yi  (J»N)  =0  YY  (J«-  :i)  -3.  *a  ( J*- 1  ,  1)  *COS?  1  *A  (  1  ,  K«-  1)  *  i  Y  (  1 ,  X  j  /  A  i  1 ,  1) 

1  >(3{J*-1,K*1)~3(J+l,1)*\(l,X'*’1)/Ail,1))'>YY(l/A<-.i> 

DYY  (J)  =:0EE1  ♦CC5F5*DYY  (J  ) 

DYY  (J  ♦■■N)  =0  YY  (J  +  ;l)  ♦COEFS 

3  0?  £3  =2.  ♦  J3rGK»  (0 .5*  (  YY  ( 1,  J)  *•  Y  Y  (  1  ,  J*-  N)  )  -DX3?/Y  Y  ( 1  ,  J+;.)  /  YY  .1  ,o  +  ;.)  ) 
IF(30?EA.J  1-  0.)  DYY  lJ*-3)  =DYY(J«-;i)  -SUPE?, 

123  CJSTIAJE 
itaTO 


Sa36DUTl:iE  PEDE3  V  (“I,  Y,  ) 

IrtPLICir  BEAE*6  (A-H,  Q-Z) 

DI0S:J3iJN  Y(8,20)  (“00)  ,AdO,20)  ,  3  (20,20) 


!¥««««  «««»«^ 


SOaPOOTIOE  CAZi.£0  3Y  DIFSOS  STOriES  THE  PARTIAL 
DESIVATI7E3  jf  THE  DIFFE3SMIAL  ECOATlOXo  PROVIDED  I.. 
S'J3ROUTI.IE  j^^c'J'A.  THE  PAPTIAi.  DERIVATIVES  (THE  JACOEIAO) 
iJEIiE  EVAlJAi  ED  10  THE  LAST  SECTlUO  10  CErtPTE?  5. 


C 

c  *  •  > 


v. 

c 

c 


1  *  **  i 


i  «  «  *  «  «  . 


COM.H  DM/ E  1 /O  OEF  1 ,  CO  EF5 
COH.IOO/E  J/.N  ,  N  T 

COtlHOki/ 1**^  /  ti,  a 

CoH.IDN/;.  10/3aI3K  ,  D  XS  ? 

Du  la  !\E  =  1  ,  0 
DO  20  il=l,.4 

P'f  (II  +  (2*  2)  *  :.•)  =CCE  J  1«CGEr  3*0  (II  <■  i  ,  XK  f  1 ) 

P  «'  ( 1 1  (  2  •  tl  K  -  1 )  *  % )  --CC  E  ?5  *  J  •  •  3  ( 1 1  ♦  1  ,  1 ;  •  C  0  E  r  1  *  X  (  1  ,  X 1 )  /  A  v  1  ,  '  i 
IE  (  II.  E^.  K  X)  DO  TC  30 
jU  .  D  4  V, 

30  3uP=2.*DjIjX*(D.5*{Y(1,IIj  ♦Y(1,ZI*X)  )  -  D.\wP/Y  ^l,II’’.i//i  ^1,  I..*-:,;  j 

I  F  (  JU  ?.  01.  0.  )  P  (  II  ♦  (  2  ■‘XK-  1)  *;. )  =P  /.  ( I  I  ^  (  2*  K  K-  i )  «  :. )  -  DEI  D  X 
(  II*2*o*H>  (2*KX-2)  *';)  =0. 


40 


1 
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♦COEr  5 

Ir(II.E^.KK)  JC  TC  uO 
uC  TO  2J 

oO  ir(  jdV.  or.  0  - )  ?w  (ii+2*M*.Sf  1)  =?„  ( ii4-2«o  (2*.\k-  i )  ».<)  -2.  * 

1  0  3iaK»  (.  b*-2.  •D^i?/T  (  1,  XI+;.)  / 

2  1 
20  CUN'TIj'JE 

10  cu.inwuE 
l\  ET  J  dO 
ENO 


THE  CTrlEii  iJBSO'ITIilES  CALLED  Si  THIS  .’SuCiSAl^  STAY  TEE  jA.lE  C 
AS  TiiJSE  ::i  TEE  PI5ST  PSOSaA.I.  C 


*««««« 


